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Here, we review the basic concepts and apphcations of the phase- field- crystal (PFC) method, 
which is one of the latest simulation methodologies in materials science for problems, where 
atomic- and microscales are tightly coupled. The PFC method operates on atomic length and 
diffusive time scales, and thus constitutes a computationally efficient alternative to molecu- 
lar simulation methods. Its intense development in materials science started fairly recently 
following the work by Elder etal. [Phys. Rev. Lett. 88 (2002), p. 245701]. Since these initial 
studies, dynamical density functional theory and thermodynamic concepts have been linked 
to the PFC approach to serve as further theoretical fundaments for the latter. In this review, 
we summarize these methodological development steps as well as the most important appli- 
cations of the PFC method with a special focus on the interaction of development steps taken 
in hard and soft matter physics, respectively. Doing so, we hope to present today's state of 
the art in PFC modelling as well as the potential, which might still arise from this method in 
physics and materials science in the nearby future. 
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Figure 1. Pattern formation on microscopic to cosmic length scales. From left to right: multiple spiralling 
nanoscale terraces starting from a central heterogeneity. (Reproduced with permission from C. Klemenz, 
Hollow cores and step-bunching effects on (001) YBCO surfaces grown by liquid-phase epitaxy, J. Cryst. 
Growth 187 (1998), 221-227, no. 2, DOT: 10.1016/S0022-0248(97)00866-X © 1998 by Elsevier.) Cellular 
slime mould self-organized into a five-arm spiral structure. (Reproduced with permission from B. Vasiev, F. 
Siegert, and C. Weijer, Multiarmed spirals in excitable media, Phys. Rev. Lett. 78 (1997), 2489-2492, no. 
12, DOT: 10.1103/PhysRevLett.78.2489 © 1997 by the American Physical Society.) Messier 100, a multi- 
arm spiral galaxy in the Virgo Supercluster, 60 million light-years from earth. (Credit: ESO/IDA/Danish 
1.5m/R. Gendler, J.-E. Ovaldsen, C. C. Thone, and C. Feron.) 

1. Introduction 

Pattern formation has been observed in complex systems from microscopic to cos- 
mic scales (for examples, see figure 1), a phenomenon that has been exciting the 
fantasy of humanity for a long time. Nonequilibrium systems in physics, chem- 
istry, biology, mathematics, cosmology, and other fields produce an amazingly rich 
and visually fascinating variety of spatiotemporal behaviour. Experiments and sim- 
ulations show that many of such systems - reacting chemicals, bacteria colonies, 
granular matter, plasmas - often display analogous dynamical behaviour. The wish 
to find the origin of the common behaviour has been driving the efforts for find- 
ing unifying schemes that allow the assigning of many of these processes into a 
few universality classes. Pattern formation and the associated nonlinear dynamics 
have received a continuous attention of the statistical physics community over the 
past decades. Reviews of the advances made in different directions are available 
in the literature and range from early works on critical dynamics [1] via phase- 
separation [2] and pattern formation in nonequilibrium systems [3, 4] to recent 
detailed treatments of the field in books (e.g., references [5-7]). 

In particular, Seul and Andelman [4] described pattern formation on the 
mesoscale as manifestation of modulated structures. Within this approach, the 
modulated phases are stabilized by competing attractive and repulsive interac- 
tions, which favour inhomogeneities characterized by a certain modulation length 
scale. The modulations are described by a single scalar order parameter. As out- 
lined in reference [4], the idea of Seul and Andelman can be applied to a large 
variety of systems ranging from Langmuir films over semiconductor surfaces and 
magnet garnets to polyelectrolyte solutions. Furthermore, the pioneering theories 
of spontaneous domain formation in magnetic materials and in the intermediate 
state of type I superconductors has been reinterpreted within this framework. 

In the past decade, special attention has been paid to a similar model, whose 
mathematical formulation has been laid down decades earlier to address hydro- 
dynamic instabilities [10] and to describe the transition to the antiferromagnetic 
state in liquid ^He or to a non-uniform state in cholesteric liquid crystals [11], 
whereas recently it has been employed for the modelling of crystallization in un- 
dercooled liquids on the atomic scale [12]. This approach is known to the materials 
science community as the phase-field-crystal (PFC) model [12], and has proven to 
be an amazingly efficient tool for addressing crystalline self-organization/pattern 
formation on the atomistic scale. 
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The PFC approach attracts attention owing to a unique situation: the crystal- 
hzation of hquids is traditionally addressed on this scale by the density functional 
theory (DFT) [13-15], whose best developed non-perturbative version, known as 
the fundamental-measure theory (FMT) [16], leads to unprecedented accuracy for 
such properties as the liquid-solid interfacial free energy [17, 18] or the nucleation 
barrier [17]. However, handling of large systems is hampered by the complexity 
of such models. In turn, the PFC model, being a simplistic DFT itself, incor- 
porates most of the essential physics required to handle freezing: it is atomistic, 
anisotropics and elasticity are automatically there, the system may choose from 
a variety of periodic states [such as body-centred cubic (bcc), face-centred cubic 
(fee), and hexagonal close packed (hep)] besides the homogeneous fluid, etc. The 
free-energy functional is fairly simple having the well-known Swift-Hohenberg (SH) 
form 



where ijj is the reduced particle density and /? a reduced temperature, while is 
the absolute value of the wave number vector the system prefers. (In equation (1) 
all quantities are dimensionless.) This together with the assumption of overdamped 
conservative (diffusive) dynamics (a major deviation from the non-conservative dy- 
namics of the SH model) leads to a relatively simple equation of motion (EOM) 
that, in turn, allows the handling of a few times 10^ atoms on the diffusive time 
scale. Such abilities can be further amplified by the amplitude equation versions 
[19] obtained by renormalisation group theory, which combined with advanced 
numerics [20] allows for the handling of relatively big chunks of material, while re- 
taining all the atomic scale physics. Such a coarse-grained PFC model, relying on 
equations of motion for the amplitudes and phases, can be viewed as a physically 
motivated continuum model akin to the highly successful and popular phase-field 
(PF) models [21-25], which however usually contain ad hoc assumptions. Accord- 
ingly, the combination of the PFC model with coarse graining establishes a link 
between DFT and conventional PF models, offering a way for deriving the latter 
on physical grounds. 

In its simplest formulation, defined above, the PFC model consists of only a sin- 
gle model parameter /3 (provided that length is measured in k^^ units). Still it has 
a fairly complex phase diagram in three spatial dimensions (3D), which has stabil- 
ity domains for the bcc, fee, and hep structures, as opposed to the single triangular 
crystal structure appearing in 2D. Introducing additional model parameters, re- 
cent extensions of the PFC model either aim at further controlling of the predicted 
crystal structure or attempt to refine the description of real materials. Other ex- 
tensions address binary systems, yet others modify the dynamics via considering 
further modes of density relaxation besides the diffusive one, while adopting a free 
energy that ensures particle conservation and allows assigning inertia to the par- 
ticles. In a few cases, PFC models tailored to specific apphcations have reached 
the level of being quantitative. Via the PFC models, a broad range of exciting 
phenomena became accessible for atomistic simulations (see Tab. 1), a situation 
that motivates our review of the present status of PFC modelling. 

While Tab. 1 contains a fairly impressive list, it is expected to be only the begin- 
ning of the model's employment in materials science and engineering. For example, 
true knowledge-based tailoring of materials via predictive PFC calculations is yet 
an open vision, for which a number of difficulties need to be overcome. We are going 
to review a few of the most fundamental ones of these open issues. For example. 




(1) 
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the PFC models still have to establish themselves as widely accepted simulation 
tools in materials engineering/ design, which requires methodological advances in 
various directions such as (a) ensuring the quantitativeness of PFC predictions for 
practically relevant (multi-phase multi-component) materials and (b) a consistent 
extension of PFC modelling to some essential circumstances such as non-isothermal 
problems, coupling to hydrodynamics, or handling of non-spherical molecules. 

So far, only limited reviews of PFC modelling are available [25, 62]. Therefore, we 
give a comprehensive overview of PFC modelling in the present review. Especially, 
we present the main achievements of PFC modelling and demonstrate the potential 
these models offer for addressing problems in physics and materials science. We pay 
special attention to the similarities and differences of development steps taken in 
hard and soft matter physics, respectively. The rest of our review article is struc- 
tured as follows: in section 2, we present a detailed theoretical derivation of the 
PFC model on the basis of dynamical density functional theory (DDFT). Section 
3 is devoted to some essential features of the PFC model and its generalisations 
including the realization of different crystal lattices, the predicted phase diagrams, 
anisotropy, and some specific issues such as glass formation, application to foams, 
and the possibility for coupling to hydrodynamics. Section 4 addresses nucleation 
and pattern formation in metallic alloys, whereas section 5 deals with the appli- 
cation of the PFC models to prominent soft matter systems. Finally, in section 6, 
we offer a few concluding remarks and an outlook to probable developments in the 
near future. 



2. From density functional theory to phase-field-crystal models 

Freezing and crystallization phenomena are described best on the most fundamen- 
tal level of individual particles, which involves the microscopic size and interaction 

Table 1. A non-exclusive collection of phenomena addressed using PFC techniques. 



Phenomena 


References 


Liquid-solid transition: 




- dendrites 


[26-30] 


- eutectics 


[26, 28, 29, 31] 


- homogeneous nucleation 


[28, 30-33] 


- heterogeneous nucleation 


[30, 31, 34, 35] 


- grain-boundary melting 


[36, 37] 


- fractal growth 


[38, 39] 


- crystal anisotropy 


[33, 38, 40-44] 


- density /solute trapping 


[38, 39, 45] 


- glass formation 


[35, 46, 47] 


- surface alloying 


[48, 49] 


Colloid patterning 


[33] 


Grain-boundary dynamics 


[50] 


Crack propagation 


[50] 


Elasticity, plasticity, dislocation dynamics 


[12, 50-54] 


Kirkendall effect 


[55] 


Vacancy transport 


[56] 


Liquid phase separation with colloid 




accumulation at phase boundaries 


[57] 


Liquid crystals 


[58-60] 


Formation of foams 


[61] 
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Figure 2. Levels of description with the corresponding methods and theories (schematic). 

length scale of the particles (see figure 2). The individual dynamics of the parti- 
cles happens correspondingly on a microscopic time scale. In the following, two 
different classes of materials, namely molecular and colloidal materials, need clear 
distinction. The former comprise metals as well as molecular insulators and semi- 
conductors. We consider these molecular systems as classical particles, where the 
quantum-mechanical nature of the electrons merely enters via effective molecular 
force fields. The corresponding molecular dynamics is governed by Newton's sec- 
ond law. Hence the length scale is atomic (about a few Angstroms) and the typical 
time scale is roughly a picosecond. 

The latter material class of colloidal systems involves typically mesoscopic par- 
ticles immersed in a molecular viscous fluid as a solvent that are interacting via 
effective forces [63] . These colloidal suspensions have a dimension typically in the 
range between a nanometre and a micrometer and are therefore classical parti- 
cles. Thus, the corresponding "microscopic" length scale describing their extension 
and interaction range is much bigger than for the molecular materials. The indi- 
vidual particle dynamics is Brownian motion [64, 65], i.e., it is completely over- 
damped^ superimposed with stochastic kicks of the solvent. The corresponding 
coarse-grained Brownian time scale upon which individual particle motion occurs 
is much longer (about a microsecond) [68]. 



-"^It is interesting to note that there are also mesoscopic particle systems with Newtonian dynamics, which 
are virtually undamped. These are realized in so-called complex plasmas [66, 67], where dust particles are 
dispersed and levitated in a plasma. 
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In terms of static equilibrium properties (such as structural correlations and 
phase transitions), both metals and colloids can just be regarded as classical in- 
teracting many-body systems. For this purpose, density functional theory (DFT) 
was developed [14, 15, 69]: DFT is a microscopic theory, i.e., it starts with the 
(effective) interparticle interactions and predicts the free energy and the static 
many-body correlations. In principle, DFT is exact, but for practical applications 
one has to rely on approximations. 

In the past years, it has become clear that DFT is an ideal theoretical framework 
to justify and to derive the free-energy functional of coarse-grained models as the 
PFC approach [26, 70]. PFC models keep the microscopic length scale, but describe 
the microscopically structured density field in a very rough way, for example, by 
keeping only its first Fourier modes for a crystal. Although some microscopic details 
are lost, the basic picture of the crystal is kept and much larger system sizes can 
be explored numerically. The PFC models are superior to simple PF models, which 
work with a single order parameter on a more coarse-grained regime. Finally, there 
are also phenomenological hydrodynamic approaches that are operating on the 
macroscopic length and time scale. 

This pretty transparent hierarchy of length scales for static equilibrium proper- 
ties gets more complex for the dynamics. In order to discuss this in more detail, 
it is advantageous to start with the colloidal systems first. Here, the individual 
dynamics is already dissipative and overdamped: the "microscopic" equations gov- 
erning the colloidal Brownian dynamics are either the Langevin equation for the 
individual particle trajectories or the Smoluchowski equation for the time evolution 
of the many-body probability density [71, 72]. Both approaches are stochastically 
equivalent [73]. In the end of the past century, it has been shown that there is 
a dynamic generalisation of DFT, the DDFT, which describes the time evolution 
of the many-body system within the time-dependent one-body density as a gen- 
eralized deterministic diffusion equation. This provides a significant simplification 
of the many-body problem. Unfortunately, DDFT is not on the same level as the 
Smoluchowski or Langevin picture since an additional adiabaticity approximation 
is needed to derive it. This approximation implies, that the one-body density is 
a slowly relaxing variable and all higher density correlations relax much faster 
to thermodynamic equilibrium [74]. Fortunately, the adiabaticity approximation 
is reasonable for many practical applications except for situations, where fluctua- 
tions play a significant role. Now, DDFT can be used as an (approximate) starting 
point to derive the dynamics of a PFC model systematically [70] . This also points 
to alternative dynamical equations, which can be implemented within a numeri- 
cally similar effort as compared to ordinary PFC equations, but are a bit closer to 
DDFT. 

For Newtonian dynamics, on the other hand, intense research is going on to derive 
a similar kind of DDFT [75-78]. Still the diffusive (or model B) dynamics for a 
conserved order-parameter field can be used as an effective dynamics on mesoscopic 
time scales with an effective friction. Then, the long-time self-diffusion coefficient 
sets the time scale of this process. One should, however, point out that the PFC 
dynamics for molecular systems is dynamically more coarse-grained than for their 
colloidal counterparts. 



2.1. Density functional theory 

Density functional theory (DFT) is a microscopic theory for inhomogeneous com- 
plex fluids in equilibrium [14, 15, 69, 79] that needs only the particle interactions 
and the underlying thermodynamic conditions as an input. The central idea is to 
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express the free energy of the many-body system as a functional of the inhomoge- 
neous one-body density. As it stands originally, DFT is a theory for static quanti- 
ties. Most of the actual applications of DFT are for spherically symmetric pairwise 
interactions between classical particles (mostly hard spheres) [14, 69, 80], but they 
can also be generalized to anisotropic interactions (as relevant for non-spherical 
hard bodies or molecules) [81-84]. One of the key applications of DFT concerns 
the equilibrium freezing and melting [14, 15, 85, 86] including the fluid-solid inter- 
face [87-90]. Further information about DFT and a detailed historic overview can 
be found in several articles and books like references [87, 91-94]. 

More recently, static DFT was generalized towards time-dependent processes 
in nonequilibrium. The extended approach is called dynamical density functional 
theory (DDFT). DDFT was first derived in 1999 for isotropic Brownian particles 
by Marconi and Tarazona [95, 96] starting from the Langevin picture of individual 
particle trajectories. An alternate derivation based on the Smoluchowski picture 
was presented in 2004 by Archer and Evans [97]. In both schemes an additional 
adiabaticity approximation is needed: correlations of high order in nonequilibrium 
are approximated by those in equilibrium for the same one-body density. These 
derivations were complemented by a further approach on the basis of a projection 
operator technique [74]. The latter approach sheds some light on the adiabaticity 
approximation: it can be viewed by the assumption that the one-body density 
relaxes much slower than any other density correlations of higher order. DDFT 
can be flexibly generalized towards more complex situations including mixtures 
[98], active particles [99], hydrodynamic interactions [100, 101], shear flow [102] 
and non-spherical particles [103, 104]. However, as already stated above, it is much 
more difficult to derive a DDFT for Newtonian dynamics, where inertia and flow 
effects invoke a treatment of the momentum density field of the particles [75-78]. 

In detail, DFT gives access to the free energy for a system of N classical par- 
ticles, whose centre-of-mass positions are defined through the vectors with 
i G {1,...,A^}, by the one-particle density p(r), which provides the probability 
to find a particle at position r. Its microscopic definition is 



with the normalised classical canonical (or grand canonical) ensemble- average ( • ). 
At given temperature T and chemical potential /i, the particles are interacting via 
a pairwise (two-body) potential U2{ri — r2). Furthermore, the system is exposed to 
an external (one-body) potential Ui{r) (describing, for example, gravity or system 
boundaries), which gives rise in general to an inhomogeneous one-particle density 
p(r). DFT is based on the following variational theorem: 

There exists a unique grand canonical free- energy functional ^{T, Jjl, [p{'^)]) of the one- 
particle density p{r), which becomes minimal for the equilibrium one-particle density 



If the grand canonical functional Vt{T^p^ [p{^)]) i^ evaluated at the equilibrium one- 
particle density p{r), it is the real equilibrium grand canonical free energy of the 
inhomogeneous system. 

Hence, DFT establishes a basis for the determination of the equilibrium one-particle 
density field p(r) of an arbitrary classical many-body system. However, in practice, 




i=l 



(2) 



p(r).- 



5n{T,^i, [p{v)]) 
5p{r) 



= 0. 



(3) 
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the exact form of the grand canonical free-energy density functional 0(r, /i, [p(r)]) 
is not known and one has to rely on approximations. Via a Legendre transform, the 
grand canonical functional can be expressed by an equivalent Helmholtz free-energy 
functional J^{T, [p(r)]), 

0(r,/i, [p(r)]) = J-(r, [p(r)]) - /i ydrp(r) , (4) 

with V denoting the system volume. The latter is conveniently split into three 
contributions: 

T{T, [p(r)]) = J-id(r, [p(r)]) + T,.,{T, [p(r)]) + T,.t{T, [p(r)]) . (5) 

Here, ^id(r, [/>(r)]) is the (exact) ideal gas free- energy functional [69] 

J-id(T, [p(r)]) = fceT Jdr p{r) {ln{A^ p{r)) - l) , (6) 

where /cb is the Boltzmann constant and A the thermal de Broglie wavelength. 
The second term on the right-hand-side of equation (5) is the excess free-energy 
functional Te^dT, [p(r)]) describing the excess free energy over the exactly known 
ideal-gas functional. It incorporates all correlations due to the pair interactions 
between the particles. In general, it is not known explicitly and therefore needs 
to be approximated appropriately [14, 69]. The last contribution is the external 
free- energy functional [69] 

J-e,t(T,[p(r)]) = Jdrp{r)Uiir). (7) 

A formally exact expression for J^exc(7'Jp(r)]) is gained by a functional Taylor 
expansion in the density variations Ap(r) = p(r) — pref around a homogeneous 
reference density pref [69, 85]: 

CXD ^ 

J^UT, [p{r)]) = miPrei) + k^TY,-,:Ftl{T, (8) 



n=l 



with 



J-i^)(T, [p(r)]) = - /dn • • • /dr„ c(-)(ri, . . . , r„) J] Ap(r,) . (9) 

i=l 

Here, c(^)(ri, . . . ,r^) denotes the nth-order direct correlation function [92] in the 
homogeneous reference state given by 



cW(ri,...,r„) = 



(10) 

Pref 



depending parametrically on T and pref • 
In the functional Taylor expansion (8), the constant zeroth-order contribution 
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is irrelevant and the first-order contribution corresponding to n = 1 is zero.^ The 
higher-order terms are nonlocal and do not vanish in general. 

In the simplest nontrivial approach, the functional Taylor expansion is trun- 
cated at second order. The resulting approximation is known as the Ramakrishnan- 
Yussouff theory 



:F,^,{T, [p{v)\) = -^/cbT jdvi jdv2 c(2)(ri - r2)Ap(ri)Ap(r2) (11) 

and predicts the freezing transition of hard spheres both in three [85] and two 
spatial dimensions (2D) [105].^ The Ramakrishnan-Yussouff approximation needs 
the fluid direct pair-correlation function c^'^\vi — Y2) as an input. For exam- 
ple, c(^)(ri — r2) can be gained from liquid integral equation theory, which links 
c^'^)[yi — Y2) to the pair-interaction potential C/2(ri — Y2). Well-known analytic ap- 
proximations for the direct pair-correlation function include the second-order virial 
expression [109] 

.(2)/ A U2{yi-Y2] 



c^'\v, - r,) = exp(^- ^YbT ) ~ ^ ' ^^^^ 

The resulting 0ns ager functional for the excess free energy becomes asymptotically 
exact in the low density limit [92]. An alternative is the random-phase or mean-field 
approximation 

(2)/ N f^2(ri-r2) 
\ri - rs) = — . (13) 

For bounded potentials, this mean-field approximation becomes asymptotically ex- 
act at high densities [103, 110-112]. Non-perturbative expressions for the excess 
free-energy functional for colloidal particles are given by weighted- density approx- 
imations [81, 84, 87, 113, 114] or follow from FMT [82, 115]. FMT was originally 
introduced in 1989 by Rosenfeld for isotropic particles [16, 86, 91, 116] and then 
refined later [80, 117] - for a review, see reference [79]. For hard spheres, FMT 
provides an excellent approximation for the excess free-energy functional with 
an unprecedented accuracy. It was also generalized to arbitrarily shaped particles 
[82, 115, 118]. 



2.2. Dynamical density functional theory 

2.2.1. Basic equations 

Dynamical density functional theory (DDFT) is the time-dependent analogue 
of static DFT and can be classified as linear-response theory. In its basic form, 
it describes the slow dissipative nonequilibrium relaxation dynamics of a system 
of N Brownian particles close to thermodynamic equilibrium or the behaviour in 
a time-dependent external potential ?7i(r,t). Now a time-dependent one-particle 



^This follows from the representation (9) under consideration of the translational and rotational symmetries 
of the isotropic bulk fluid that also apply to the direct correlation function c(^)(ri) = const. 
^More refined approaches include also the third-order term [106] with an approximate triplet direct corre- 
lation function [107, 108]. 
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density field is defined via 

N 



p(r,t) = (5]%-r,(i))), (14) 



i=l 

where ( • ) denotes the normalised classical canonical noise-average over the particle 
trajectories and t is the time variable. 

This one-particle density is conserved and its dynamics is assumed to be dissi- 
pative via the generalized (deterministic) diffusion equation 



dt k^T ' ' ' ^ 5p(v,t) 

Here, Dt denotes a (short-time) translational diffusion coefficient for the Brownian 
system. Referring to equations (3) and (4), the functional derivative in the DDFT 
equation can be interpreted as an inhomogeneous chemical potential 

such that the DDFT equation (15) corresponds to a generalized Pick's law of 
particle diffusion. As already mentioned, DDFT was originally invented [95, 97] for 
colloidal particles, which exhibit Brownian motion, but is less justified for metals 
and atomic systems whose dynamics are ballistic [75, 76]. 

2.2.2. Brownian dynamics: Langevin and Smoluchowski picture 

The DDFT equation (15) can be derived [95] from Langevin equations that de- 
scribe the stochastic motion of the N isotropic colloidal particles in an incompress- 
ible liquid of viscosity rf at low Reynolds number (Stokes limit). In the absence of 
hydrodynamic interactions between the Brownian particles, these coupled Langevin 
equations for the positions r^(t) of the colloidal spheres with radius describe 
completely overdamped motion plus stochastic noise [71, 73]: 

ri = r'(Fi + fi) , i = l,...,N. (17) 

Here, ^ is the Stokesian friction coefficient (^ = 67rr]Rs for spheres of radius Rq 
with stick boundary conditions) and 

Fi{t) = -VrV{Ti,...,rN,t) (18) 

are the deterministic forces caused by the total potential 

C/(ri, . . . ,r7v,t) = C/ext(ri, . . . ,r7v,t) + C/int(ri, . . . , r^v) (19) 

with 



N 

[7ext (r 1 , . . . , r AT , t) = ^ /7i (r„ t) (20) 

i=l 
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and 

N 



C/int(ri,...,r^) = ^2(r^-r,) . (21) 

Z,J = 1 

i<j 

On top of these deterministic forces, also stochastic forces f^(t) due to thermal 
fluctuations act on the Brownian particles. These random forces are modelled by 
Gaussian white noises with vanishing mean values 

(f,(i))=0 (22) 

and with Markovian second moments 

(f,(ti) ® f,(t2)) = 2ikBT15ij5{ti - t2) , (23) 

where ® is the ordinary (dyadic) tensor product (to make the notation compact) 
and 1 denotes the 3 x 3-dimensional unit matrix. This modelling of the stochastic 
forces is dictated by the fluctuation-dissipation theorem, which for spheres yields 
the Stokes-Einstein relation = k^T/^^ [64], that couples the short-time diffusion 
coefficient of the colloidal particles to the Stokes friction coefficient ^. 

An alternate description of Brownian dynamics is provided by the Smoluchowski 
picture, which is stochastically equivalent to the Langevin picture [72, 73]. The 
central quantity in the Smoluchowski picture is the A/'-particle probability density 
P(ri, . . . ^Yjsf^t) whose time evolution is described by the Smoluchowski equation 
[73, 119] 

^P(ri , . . . , r^, t) = £ P(ri , . . . , r^, t) (24) 
with the Smoluchowski operator 

£ = gTEv. (v„ ^<^'- ;'-"-" +V,V (25) 



While the A^-particle probability density P(ri, . . . , r^y, t) in this Smoluchowski 
equation is a highly nontrivial function for interacting particles, it is often sufficient 
to consider one-body or two-body densities. The one-particle probability density 
P(r,t) is proportional to the one-particle number density p(r,t). In general, all 
n-particle densities with n ^ N can be obtained from the A^-particle probability 
density P(ri, . . . , r^v, ^) by integration over the remaining degrees of freedom: 

(n , . . . , r,, t) = jj/^^ Jdrn+1 • • • JdTN P(ri , . . . , r^, t) . (26) 

2. 2. 3. Derivation of DDFT 

We now sketch how to derive the DDFT equation (15) from the Smoluchowski 
picture following the idea of Archer and Evans [97]. Integrating the Smoluchowski 
equation (24) over the positions of iV — 1 particles yields the exact equation 

p{r,t) = DTVr- (Vrp(r,t) - + P^VMr,t)\ (27) 
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for the one-particle density p(r,t), where 

F{v, t) = -jdv'p^^^ (r, r', t)VrU2{v - v') (28) 

is an average force, that in turn depends on the nonequiUbrium two-particle den- 
sity />(^)(ri, r2, t). This quantity is approximated by an equilibrium expression. To 
derive this expression, we consider first the equilibrium state of equation (27). This 
leads to 

F(r) = kBTVrp{v) + p(r)VrC7i(r) , (29) 

which is the first equation of the Yvon-Born-Green hierarchy, with a "substitute" 
external potential Ui{y). In equilibrium, DFT implies 

^ ^i^(r,/x,[p(r)]) ^ snTM^)]) _ „ 

= ^BTln(A3,(r)) + ^^^^^1^ +I7,W 

and after application of the gradient operator 

^ VrP(r) , ^ ^J-exc(r,[p(r)]) , ^ ^ . . 
= fcB-T — r^ + Vr |-Vrt/i(r). (31) 

A comparison of equations (29) and (31) yields 

m - -p(r)Vr^^4^;^ • (^2) 

It is postulated, that this relation also holds in nonequilibrium. The nonequilibrium 
correlations are thus approximated by equilibrium ones at the same p(r) via a suit- 
able "substitute" equilibrium potential C/i(r). With this adiabatic approximation^ 
equation (27) becomes 

p(r, t) = i^TVr- (VrPir^t) + + ^Vrf/i(r, t)j , (33) 

which is the DDFT equation (15). 

It is important to note that the DDFT equation (15) is a deterministic equation, 
i.e., there are no additional noise terms. If noise is added, there would be double 
counted fluctuations in the equilibrium limit of equation (15) since J^(T, [p]) is 
the exact equilibrium functional, which in principle includes all fluctuations. The 
drawback of the adiabatic approximation, on the other hand, is that a system is 
trapped for ever in a metastable state. This unphysical behaviour can be changed 
by adding noise on a phenomenological level though violating the caveat noted 
above. A pragmatic recipe is to add noise only when fluctuations are needed to 
push the system out of a metastable state or to regard a fluctuating density field as 
an initial density profile for subsequent deterministic time evolution via DDFT. In 
conclusion, the drawback of the adiabatic approximation is that DDFT is some kind 
of mean-field theory. For example, DDFT as such is unable to predict nucleation 
rates. It is rather a realistic theory, if a systematic drive pushes the system, as 
occurs, for example, for crystal growth. 
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Figure 3. Crystallization starting at a colloidal cluster. The plots show DDFT results for the time- 
dependent density field. Apref = 0-7 (left column) and Apref = 0-6 (right column) at times t/rB = 
0,0.001,0.1,1 (from top to bottom) with the area A of a unit cell of the imposed crystalline seed, the 
Brownian time tb, and the lattice constant a = (2/(A/3pref ))"'^^^- For Apref = 0-7, the cluster is compressed 
in comparison to the stable bulk crystal, but there is still crystal growth possible, while the compression 
is too high for Ap^ef = 0-6. The initial nucleus melts in this case. (Reproduced from S. van Teeffelen, C. 
N. Likos, and H. Lowen, Colloidal crystal growth at externally imposed nucleation clusters, Phys. Rev. 
Lett. 100 (2008), 108302, no. 10, DOT: 10. 1103/PhysRevLett. 100. 108302 © 2008 by the American Physical 
Society.) 



2.2.4' Application of DDFT to colloidal crystal growth 

An important application of DDFT is the description of colloidal crystal growth. 
In reference [120], DDFT was applied to two-dimensional dipoles, whose dipole 
moments are perpendicular to a confining plane. These dipoles interact with a re- 
pulsive inverse power-law potential U2{'r) = UQr~^, where r = |r| is the interparticle 
distance. This model can be realized, for example, by superparamagnetic colloids 
at a water-surface in an external magnetic field [121]. Figures 3 and 4 show DDFT 
results from reference [120]. 

In figure 3, the time evolution of the one-particle density of an initial colloidal 
cluster of 19 particles arranged in a hexagonal lattice is shown. This prescribed 
cluster is surrounded by an undercooled fluid and can act as a nucleation seed, if 
its lattice constant is chosen appropriately. The initial cluster either initiates crystal 
growth (left column in figure 3) or the system relaxes back to the undercooled fluid 
(right column). 

A similar investigation is also possible for other initial configurations like rows of 
seed particles. Figure 4 shows the crystallization process starting with six infinitely 
long particle rows of a hexagonal crystal, where a gap separates the first three rows 
from the remaining three rows. If this gap is not too big, the density peaks rearrange 
and a growing crystal front emerges. 
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Figure 4. Crystallization starting at two triple-rows of hexagonally crystalline particles that are separated 
by a gap. The contour plots show the density field of a growing crystal at times I/tb = 0, 0.01, 0.1, 0.63, 1 
(from top to bottom). (Reproduced from S. van Teeffelen, C. N. Likos, and H. Lowen, Colloidal crystal 
growth at externally imposed nucleation clusters^ Phys. Rev. Lett. 100 (2008), 108302, no. 10, DOI: 
10. 1103/PhysRevLett. 100. 108302 © 2008 by the American Physical Society.) 

2.3. Derivation of the PFC model for isotropic particles from DFT 

Though approximate in practice, DFT and DDFT can be regarded to be a high 
level of microscopic description, which provides a framework to calibrate the more 
coarse-grained PFC approach. In this section, we at first describe the derivation for 
spherical interactions in detail and then focus more on anisotropic particles. There 
are two different aspects of the PFC modelling, which can be justified from DFT 
respectively DDFT, namely statics and dynamics. The static free energy used in 
the PFC model was first derived from DFT by Elder and co-workers in 2007 [26], 
while the corresponding dynamics was derived from DDFT by van Teeffelen and 
co-workers in 2009 [70]. We follow the basic ideas of these works in the sequel. 

2.3.1. Free-energy functional 

For the static part, we first of all define a scalar dimensionless order-parameter 
field ipijcY by the relative density deviation 

p(r) = pref(l + ^(r)) (34) 

around the prescribed fluid reference density pref- This relative density deviation 
'^(r) is considered to be small, |'0(r)| <C 1, and slowly varying in space (on the 
microscale). The basic steps to derive the PFC free energy are threefold: i) insert the 
parametrisation (34) into the (microscopic) free-energy functional (5), ii) Taylor- 
expand systematically in terms of powers of '^(r), iii) perform a gradient expansion 
[69, 88, 122-124] of '^(r). Consistent with the assumption that density deviations 
are small, the Ramakrishnan-Yussouff approximation (11) is used as a convenient 
approximation for the free-energy functional. 



■"^ Notice that the order-parameter field '0(r) introduced here is not identical with the field il^ir) in equation 
(1), although both fields are dimensionless. 
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For the local ideal gas free-energy functional (6) this yields^ 

J-idlV'lr)] = Fo + pref k^T (^V + Y - Y + ^) (35) 

with the irrelevant constant Fq = p^Q^V k^T {hi{h^ p^^f) — 1). The Taylor expansion 
is performed up to fourth order, since this is the lowest order which enables the for- 
mation of stable crystalline phases. The nonlocal Ramakrishnan-Yussouff approxi- 
mation (11) for the approximation of the excess free-energy functional ^exc['^(i")] is 
gradient-expanded to make it local. For this purpose, it is important to note that 
- in the fluid bulk reference state - the direct pair-correlation function c^'^\yi —Y2) 
entering into the Ramakrishnan-Yussouff theory has the same symmetry as the 
interparticle interaction potential U2{ti —^2). For radially symmetric interactions 
(i.e., spherical particles), there is both translational and rotational invariance im- 
plying 

c(^) (n , r2) = c(^) (ri -V2) = c(^) (r) (36) 

with the relative distance r = |ri — r2|. Then, as a consequence of equation (36), the 
Ramakrishnan-Yussouff approximation is a convolution integral. Consequently, a 
Taylor expansion of the Fourier transform c^^^ (k) of the direct correlation function 
in Fourier space (around the wave vector k = 0) 

c(2)(k)=42) + 4V + clV + ... (37) 

(2) 

with expansion coefficients q ^ becomes a gradient expansion in real space 

c(^)(r) = 4^)-cfv? + cfv^T--- (38) 

(2) 

with the gradient expansion coefficients c\ ' . Clearly, gradients of odd order vanish 
due to parity inversion symmetry c(^)(— r) = c(^)(r) of the direct pair-correlation 
function. 

The gradient expansion up to fourth order is the lowest one that makes stable 
periodic density fields possible. We finally obtain 

Te^cmr)] = Fexc - ^kBT Jdv [A^^jj^ + A^^jj^l^jj + AsV^V^V^) (39) 
with the irrelevant constant Fgxc = *^exc(Pref) and the coefficients 

roo o 1*00 roo 

A, = 47rp,ef /drr2c(2)(r) , A2 = ^Trp^ef krr'c^^\r) , A3 = ^ drr'c^'^r) 
Jo Jo Jo 

(40) 

that are moments of the fluid direct correlation function c^'^\r). 
Finally, the external free-energy functional (7) can be written as 



-^ext bPir)] = Fext + Pref jdvi,{v)Ui (r) (41) 



^This Taylor approximation of the logarithm has the serious consequence that the non-negative-density 
constraint p(r) ^ gets lost in the PFC model. 
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with the irrelevant constant Fgxt = Pref /dr C/i(r). We add as a comment here that 
this external part is typically neglected in most of the PFC calculations. Altogether, 
we obtain 

Tmr)] = pref kBT Jdr (^^iV' + ^2V'V?V + ^a^V^V' - ^ + I^) ^"^^^ 

for the total Helmholtz free-energy functional and the scaled coefficients 

A[ = ^{l-A,), A = -^A2, A = -^As (43) 

are used for abbreviation, where the coefficient A2 should be positive in order to 
favour non-uniform phases and the last coefficient is assumed to be positive 
for stability reasons. By comparison of equation (42) with the original PFC model 
(1), that was initially proposed on the basis of general symmetry considerations in 
reference [12], analytic expressions can be assigned to the unknown coefficients in 
the original PFC model: when we write the order-parameter field in equation (1) 
as '0(r) = q:(1 — 2'0(r)) with a constant a and neglect constant contributions as 
well as terms linear in ip{r) in the free-energy density, we obtain the relations 




between the coefficients in equations (1) and (42). 

2.3.2. Dynamical equations 

We turn to the dynamics of the PFC model and derive it here from DDFT. 
Inserting the representation (34) for the one-particle density field into the DDFT 
equation (15), we obtain for the dynamical evolution of the order-parameter field 



at 



= L>T Vr • ({l + V) Vr (2^; V + ^A'^Vl^P + 2A'^Vt7p - ^ + . (45) 



This dynamical equation (called PFCl model in reference [70]) still differs from 
the original dynamical equation of the PFC model. The latter can be gained by 
a further constant-mobility approximation (CMA), where the space- and time- 
dependent mobility DTP(r,t) in the DDFT equation is replaced by the constant 
mobility DxPref • The resulting dynamical equation (called PFC2 model in reference 
[70]) coincides with the original PFC dynamics given by 



at 



= D^Vl [2^^^ + 2A',Vl,p + 2A',Vt,p - Y + y) (46) 



for the time-dependent translational density ip{r, t). We remark that this dynamical 
equation can also be derived from an equivalent dissipation functional 9^ known 
from linear irreversible thermodynamics [125-127]. A further transformation of this 
equation to the standard form of the dynamic PFC model will be established in 
section 3.1.1. 
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Figure 5. Colloidal crystal growth within DDFT (upper panel) and the PFCl model (lower panel). The 
crystallization starts with an initial nucleus of 5 and 11 rows of hexagonally crystalline particles, respec- 
tively. The density field of the growing crystal is shown at times t/rB = 0,0.5, 1, 1.5. (Reproduced from 
S. van Teeffelen, R. Backofen, A. Voigt, and H. Lowen, Derivation of the phase- field- crystal model for 
colloidal solidification, Phys. Rev. E 79 (2009), 051404, no. 5, DOT: 10. 1103/PhysRevE.79. 051404 © 2009 
by the American Physical Society.) 

2.3.3. Colloidal crystal growth: DDFT versus PFC modelling 

Results of the PFCl model, the PFC2 model, and DDFT are compared for 
colloidal crystal growth in reference [70]. Figures 5-7 show the differences for the 
example of a growing crystal front starting at the edge of a prescribed hexagonal 
crystal. The underlying colloidal systems are the same as in section 2.2.4. In figure 
5, the time evolution of the one-particle density is shown for DDFT and for the 
PFCl model. The PFC2 model leads to results very similar to those for the PFCl 
model and is therefore not included in this figure. 

Two main differences in the results of DDFT and of the PFCl model are obvious: 
first, the density peaks are much higher and narrower in the DDFT results than 
for the PFCl model. While these peaks can be approximated by Gaussians in the 
case of DDFT, they are much broader sinusoidal modulations for the PFCl model. 
Secondly, also the width of the crystal front obtained within DDFT is considerably 
smaller than in the PFC approach. 

These qualitative differences can also be observed in figure 6. There, the laterally 
averaged density px^i^i^t) = {p{r^t))x2 associated with the plots in figure 5 is 
shown, where {•)x2 denotes an average with respect to X2. A further comparison 
of DDFT and the PFC approaches is possible with respect to the velocity Vf of the 
crystallization front. The corresponding results are shown in figure 7 in dependence 
of the total coupling constant F = uqv^/'^/ {k^T) and the relative coupling constant 
AF = F — Ff, where Ff denotes the coupling constant of freezing. Due to the power- 
law potential of the considered colloidal particles, their behaviour is completely 
characterized by the dimensionless coupling parameter F. When plotted versus 
AF, the growth velocity of the PFCl model is in slightly better agreement than 
that of the PFC2 model. 



3. Phase-field-crystal modelling in condensed matter physics 

The original PFC model has the advantage over most other microscopic techniques, 
such as molecular-dynamics (MD) simulations, that the time evolution of the sys- 
tem can be studied on the diffusive time scale making the long-time behaviour and 
the large-scale structures accessible [12, 50]. As already outlined in section 2.2, 
we note that the diffusion-controlled relaxation dynamics the original PFC model 
assumes is relevant for micron-scale colloidal particles in carrier fluid [70, 120], 
where the self-diffusion of the particles is expected to be the dominant way of den- 
sity relaxation. For normal liquids, the hydrodynamic mode of density relaxation 
is expected to dominate. The modified PFC (MPFC) model introduces linearised 
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Figure 6. Comparison of DDFT (upper panel) and PFCl (lower panel) results. For an analogous situation 
as in figure 5, this plot shows the laterally averaged density {xi, t) = {p(r, t))x2 at t = tb- (Reproduced 
from S. van Teeffelen, R. Backofen, A. Voigt, and H. Lowen, Derivation of the phase- field- crystal model 
for colloidal solidification, Phys. Rev. E 79 (2009), 051404, no. 5, DOI: 10.1103/PhysRevE.79.051404 © 
2009 by the American Physical Society.) 
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Figure 7. Comparison of DDFT, the PFCl model, and the PFC2 model [70]. The plot shows the velocity 
Vf of a crystallization front in the (ll)-direction in dependence of the relative coupling constant AF = F — Ff 
with the total coupling constant F and the coupling constant of freezing Ff . In the inset, the velocity 
is shown in dependence of F. (Reproduced from S. van Teeffelen, R. Backofen, A. Voigt, and H. Lowen, 
Derivation of the phase- field- crystal model for colloidal solidification, Phys. Rev. E 79 (2009), 051404, no. 
5, DOI: 10. 1103/PhysRevE.79. 051404 © 2009 by the American Physical Society.) 

hydrodynamics, realized via incorporating a term proportional to the second time 
derivative of the particle density into the EOM [51, 128], yielding a two-time-scale 
density relaxation: a fast acoustic process in addition to the long time diffusive 
relaxation of the original PFC model. A three-time-scale extension incorporates 
phonons into the PFC model [129, 130]. Another interesting group of models have 
been obtained by coarse-graining the PFC approaches [19, 20, 131], leading to 
equations of motion that describe the spatiotemporal evolution of the Fourier am- 
plitudes and the respective phase information characterizing the particle density 
field. Combined with adaptive grid schemes, the amplitude equation models are 
expected to become a numerically especially efficient class of the PFC models of 
crystallization [20]. 

Finally, we address here recent advances in the modelling of molecules or liq- 
uid crystalline systems, which are composed of anisotropic particles. There is a 
large number of molecular and colloidal realizations of these non-spherical parti- 
cles. The simplest non-spherical shape is rotationally symmetric about a certain 
axis (like rods, platelets, and dumbbells) and is solely described by an additional 
orientation vector. Liquid crystalline systems show an intricate freezing behaviour 
in equilibrium, where mesophases occur, that can possess both orientational and 
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translational ordering. Here, we show that the microscopic DFT approach for hquid 
crystals provides an excellent starting point to derive PFC-type models for liquid 
crystals. This gives access to the phase diagram of liquid crystalline phases and 
to their dynamics promising a flourishing future to predict many fundamentally 
important processes on the microscopic level. 

3.1. The original PFC model and its generalisations 

The original PFC model has several equivalent formulations and extensions that we 
review in this section. We first address the single-component PFC models. Then, an 
overview of their binary generalisations will be given. In both cases, complementing 
section 2, we start with presenting different forms of the free-energy functional, 
followed by a summary of specific forms of the EOM and of the Euler-Lagrange 
equation (ELE). Finally, we review the numerical methods applied for solving the 
EOM and ELE as well as various approaches for the amplitude equations. 

3.1.1. Single- component PFC models 
3.1.1.1. The free energy. 

The single-mode PFC model. The earliest formulation of the single-mode PFC 
(IM-PFC) model [12, 50] has been derived as a SH model with conserved dynamics 
to incorporate mass conservation. Accordingly, the dimensionless free energy of 
the heterogeneous system is given by the usual SH expression (1). We note that in 
equations (1) and (53) the analogous quantities differ by only appropriate numerical 
factors originating from the difference in the length scales. 

As already outlined in section 2.3, the free energy of the earliest and simplest 
PFC model [12] has been re-derived [26] from that of the perturbative DFT of 
Ramakrishnan and Yussouff [85], in which the free-energy difference AT = T — T^ 
of the crystal relative to a reference liquid of particle density p^ef ^^nd free energy T 
is expanded with respect to the local density difference Ap(r) = p(v) — p^^i., while 
retaining the terms up to the two-particle term (see section 2.3.1): 



dr (^pln(^-^^ -Ap^ ydriydr2Ap(ri)c(2)(ri,r2)Ap(r2) + ... 

(47) 



Fourier expanding the particle density, one finds that for the solid ps = pref (l + ^s + 
X^K exp(iK • r)), where r]s is the fractional density change upon freezing, while 
K are reciprocal lattice vectors (RLVs) and Ak the respective Fourier amplitudes. 
Introducing the reduced number density = (p — Pref)/Pref = Vs + ^k exp(iK- 
r) one obtains 



^= /dr((l + ^)ln(l + ^)-^) 

ydr iydr2 (r 1 ) c^^^ ( | r i - r2 1 ) (r2 



Pref 

Pref 



(48) 



Expanding next c^^^dri — r2|) in Fourier space, c^'^^k) ^ Cq^^ -\- c^"^ k'^ -\- c"^"^ k^ -\ , 

where c^^\k) has its first peak at = 2ix/R^^ the signs of the coefficients alternate. 



^To keep the notation simple, we ignore T and write T instead of AJ^ throughout this article. 
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(Here, i?p is the inter-particle distance.) Introducing the dimensionless two-particle 
direct correlation function c{k) = prefc^^^(/c) ^ YjjLo^'^j^'^^ = YjjLo^2j{kRp)'^^ , 
which is related to the structure factor as S'(fc) = 1/(1 — c(fc)), and integrating the 
second term on the right-hand-side of equation (48) with respect to r2 and finally 
replacing ri by r, the free-energy difference can be rewritten as 



Pref ^B^ 



Jdr (^(1 + ^) ln(l + ^)-^-t(^ ^(_i).c2, V?') ^ 



(49) 

j=0 



The reference liquid is not necessarily the initial liquid. Thus, we have here two 
parameters to control the driving force for solidification: the initial liquid number 
density po (corresponding to a reduced initial density of ipo) and the temperature 
T, if the direct correlation function depends on temperature. Taylor-expanding 
ln(l + ip) for small ip one obtains 



Pref ^B^ 



For m = 2, corresponding to the earliest version of the PFC model [12], and taking 

the alternating sign of the expansion coefficients of cf"^ into account, equation (49) 
transforms to the following form: 



Jdr (^(1 + \bo\) + f (|62|i?^V? + \b,\Rtvt)i. - T + I^) • (51) 



Pref ^bT 

Introducing the new variables 

Si = 1 + l&ol = 1 — Co [= (l/K:)/(pref ^bT), where k is the compressibility], 

= |&2p/(4|64|) [= K/{p^^i ^bT), where K is the bulk modulus], 

R = i?p (2 1641/162!)''^^^ [= the new length scale {x = Rx)^ which is now related to 

the position of the maximum of the Taylor expanded c^^^(/c)], 

and a multiplier v for the '0^-term (that accounts for the zeroth-order contribution 
from three-particle correlations), one obtains the form used by Berry et al. [36, 46]: 



= Pref ^bT Jdr (^^(A + Bs{2R^Vl + R^Vt))^ - v 



^3 ^4 
T 12 



(52) 



Here, f[ip] denotes the full (dimensional) free-energy density. 

The Swift-Hohenberg-type dimensionless form. Introducing a new set of variables, 
X = t/j = (3Ss)-'^/^'?/S, T = (Spref kBTR^B^)T^ where d is the number of spatial 
dimensions, the free energy can be transcribed into the following dimensionless 
form: 

^ = /dr(|(-e+(l + Vi)V + pf + (53) 

Here, p = -{v/2){3/ B,)^/^ = -^;(3|64|/|&2p)^/^ and e = -AB/Bs = -((1 + 
|^o|)/(|fc2p/(4|64|)) - 1), while ip = ip/{SBs)^/'^. The quantities involved in equa- 
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tion (53) are all dimensionless. Using the appropriate length unit, this expression 
becomes equivalent to equation (1) for ^ = 0. 

Equation (53) suggests that the m = 2 PFC model contains only two dimension- 
less similarity parameters, e and composed of the original model parameters. 
We note finally that even the third-order term can be eliminated. In the respective 
y = SH model (1), the state (e' = e+p^/3, ^jj' = V^-|-p/3) corresponds to the state 
(e, '0) of the original p ^ model. This transformation leaves the grand canoni- 
cal potential difference, the ELE (see section 3.1.1.3), and the EOM (see section 
3.1.1.2) invariant. Accordingly, it is sufficient to address the case p = 0. 

We stress here that in these models the approximation for the two-particle di- 
rect correlation function leads to a well defined wavelength of the density waves 
the system tends to realize (hence the name "single- mode PFC" (IM-PFC) model). 
Accordingly, any periodic density distribution that honours this wavelength repre- 
sents a local minimum of the free energy. Indeed, the IM-PFC model has stability 
domains for the bcc, fee, and hep structures (see section 3.1.1.1). Furthermore, 
elasticity and crystal anisotropics are automatically present in the model. 

Model parameters of the SH formulation have been deduced to fit to the prop- 
erties of bcc Fe by Wu and Karma [40] . 

The two-mode PFC model. An attempt has been made to formulate a free 
energy that prefers the fee structure at small e values [132], where a linear elastic 
behaviour persists. To achieve this, two well defined wavelengths were used (first 
and second neighbour RLVs), hence the name ^^two-mode PFC^^ (2M-PFC) model. 
The respective free-energy functional contains two new parameters: 

^ = /df(|(-e+(l + Vi)^(i?i + (4i + Vif))V'+f). (54) 

Here, i?i controls the relative stability of the fee and bcc structures, while k^Q\ is 
the ratio of the two wave numbers (fcrei = 2/\/3 for fee using the (111) and (200) 
RLVs). Remarkably, the IM-PFC model can be recovered for Ri oo. Model 
parameters have been deduced to fit to the properties of fee Fe by Wu and Karma 
[132]. 

We note finally that the IM-PFC and 2M-PFC models can be cast into a form 
that interpolates between them by varying a single parameter \ = Ri/{1 + Ri) £ 
[0,1] as follows [35]: 

^ = h (I (" ' + + ^'^'^^ + " + ^'^'P + t) • (^^^ 

This expression recovers the IM-PFC model limit for A = 1 {Ri oo). 

The eighth- order fitting PFC model. To approximate real bcc materials better, 
an eighth-order expansion of the Fourier transform of the direct correlation function 
around its maximum [k = fcm) has been performed recently, leading to what is 
termed the eighth-order fitting version of the phase-field-crystal (EOF-PFC) model 
[133]: 

c(2) (k) ^ c(2) (k^) - r (f - 1) ' - (I - 1) ' . (56) 
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The expansion parameters were then fixed so that the position, height, and the 
second derivative of c^'^\k) are accurately recovered. This is ensured by 

r^_ fegx(g^y(fcm) £;B=c(2)(ifc^)-c(2)(0)-r. (57) 

8 

With this choice of the model parameters and using relevant data for Fe from 
reference [40], they reported a fair agreement with MD simulation results for the 
volume change upon melting, the bulk moduli of the liquid and solid phases, and 
for the magnitude and anisotropy of the bcc-liquid interfacial free energy [133]. 

Attempts to control the crystal structure in PFC models. Greenwood and co- 
workers [134, 135] (GRP-PFC model) have manipulated the two-particle direct 
correlation function so that its peaks prefer the desired structural correlations - an 
approach that enables them to study transitions between the bcc, fee, hep, and sc 
structures. Wu, Plapp, and Voorhees [136] have investigated the possibility to con- 
trol crystal symmetries within the PFC method via tuning nonlinear resonances. 
They have proposed a general recipe for developing free-energy functional that 
realize coexistence between the liquid and periodic phases of desired crystal sym- 
metries, and have illustrated this via presenting a free-energy functional that leads 
to square-lattice- liquid coexistence in 2D. A possible extension of the method to 
the 3D case for simple cubic (sc) structures has also been discussed. 

The vacancy PFC model. The vacancy PFC (VPFC) model is an important 
extension of the PFC model that adds a term to the free energy that penalizes the 
negative values of the particle density, allowing thus for an explicit treatment of 
vacancies [56]: 

^ = h (I (" ^ + + ^'^y + T + '^(i^'i - ^')) • 

Here, h in the last term on the right-hand-side is a constant. The new term is a 
piecewise function that is zero for '0 > and positive for '0 < 0. It is then possible 
to obtain a mixture of density peaks (particles) and vacant areas (where V; ^ 0), 
resembling thus to snapshots of liquid configurations or crystalline structures with 
defects. This allows structural modelling of the fluid phase and is an important 
step towards combining the PFC model with fluid flow. The same approach has 
been used to address the dynamics of glasses [47] . 

The anisotropic PFC model. Recently, Prieler etal. [137] have extended the 
PFC approach by replacing the Laplacian in equation (1) by more general differ- 
ential operators allowing spatial anisotropy. Doing so and setting r = — (feg ~ /^) 
one arrives at the dimensionless free-energy functional of the so-called anisotropic 
PFC (APFC) model: 

Jdr ^ 2 ^ + dij Q^.Q^^ + ^'^jki Q^^Q^.Q^^Q^^ ^ ^ + ^ ^ • (59) 

Here, a^j is a symmetric matrix and hijki is a tensor of rank 4 with the symmetry 
of an elastic tensor: i ^ j^k ^ (hj) ^ (k^l) [137]. Choudhary etal. [138, 139] 
proved that based on a functional of the form (59) further crystal lattices can be 
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assessed as hexagonal, bcc, and corresponding sheared structures, for which they 
have presented the elastic parameters and identified the stationary states. 

3.1.1.2. The equation of motion. In the PFC models, different versions of the 
EOM have been employed. In all cases, conservative dynamics is assumed on the 
ground that mass conservation needs to be satisfied. (The original SH model differs 
from the IM-PFC model only in the EOM, for which the SH model assumes non- 
conserved dynamics.) Most of the PFC models rely on an overdamped conservative 
EOM (see, e.g., references [12, 26-39, 43, 46, 50, 70, 132, 137, 140]). In fact, this 
means that the particle density relaxes diffusively, a feature more characteristic to 
colloidal systems than to molten metals: in colloidal systems of particles floating 
in a carrier fluid, Brownian motion is the dominant mechanism of particle motion, 
whose properties are captured reasonably well by overdamped dynamics [38, 70, 
120] (see also section 2.2). On the contrary, in molten metals a density deficit can 
be reduced by a hydrodynamic flow of particles. Apparently, a proper dynamics 
of the solid (presence of phonons) requires three time scales [129, 130]. While 
the overdamped model has only the long diffusive time scale, the MPFC model 
realizes linearised hydrodynamics via adding a term proportional to the second 
time derivative of the density field. This new term leads to the appearance of an 
acoustic relaxation of the density (not true acoustic phonons) on a fast time scale, 
in addition to the slow diffusive relaxation at later stages [51, 128]. In a recent 
work, phonon dynamics, that acts on a third time scale, has also been introduced 
into the PFC model [129]. It has been shown that there exists a scale window, 
in which the longitudinal part of the full three-scale model reduces to the MPFC 
model, whereas the linearised hydrodynamics of the latter converges to the diffusive 
dynamics of the original PFC model for sufficiently long times [129, 130]. 

The overdamped equation of motion. In the majority of the PFC simulations, an 
overdamped conserved dynamics is assumed [that is analogous to the DDFT EOM 
(15) for colloidal systems, however, assuming here a constant mobility coefficient, 
= pQD^/{k^T)\. Accordingly, the (dimensional) EOM has the form 

|=V,.(m,V,|)+(,. (60) 

where C,p stands for the fluctuations of the density flux, whose correlator reads as 
(Cp(r,t)Cp(r',tO) = -2MpkBTVl5{Y - Y')5{t - t'). (For a discretised form of the 
conserved noise see reference [141].) 

Changing from variable p to '0, introducing — ((1 + '^/^o)^T/(Pref ^B?")), 
scaling time and distance diS t = rt and x = i?pX, where r = Rp/ (Dt(1 + '^o))? ^.nd 
inserting the free energy from equation (51), one obtains the following dimensionless 
EOM: 




with (C^(f,t)C^(f^?)) = -(2/(prefi?^))V?5(f -f05(t-tO- Analogously, the EOM 
corresponding to equation (52) has the form 
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where ((^(r, i)C;(r', t')) = -2M^kBTVlS{v - v')5{t - t'). 

Dimensionless form in Swift- Hohenherg fashion. Introducing the variables t = 
rt, X = Rx, and ^ = (3Ss)^/^^ = (3Ss)^/^('^' -2^/3) into equation (62), where 
r = R'^/{BsM^Pyq{ k^T), the EOM can be written in the form 

1^ = Vi ((- e'+{l + Vlfy + + C , (63) 

where e' = e + p^S = -{AB - {v/2f)/B, = -((1 + |6o|)/(|fc2|V(4|&4|)) - (1 + 

i;^(|64|/|62p))) and the dimensionless noise strength is o: = 2/(35g pref^^) = 
2^~^/^\h4\^~^/^ / {'iR^p^Qi |^2|4-(i/2^^ while the correlator for the dimensionless noise 
reads as (C(f , t)C(f', ?)) = -c^V?5(f - Y')5{t-t'). 

Summarizing, the dynamical m = 2 IM-PFC model has two dimensionless simi- 
larity parameters e' and a composed of the original (physical) model parameters. 
This is the generic form of the m = 2 IM-PFC model; some other formulations 
[36, 46] can be transformed into this form. 



The modified PFC model. Acoustic relaxation has been partly incorporated 
by applying an underdamped EOM. Stefanovic, Haataja, and Provatas [51] have 
incorporated a second order time derivative into the EOM of their modified PFC 
(MPFC) model, which extends the previous PFC formalism by generating dynamics 
on two time scales: 

where k and A are constants. At early times, molecular positions relax fast, consis- 
tently with elasticity theory, whereas at late times diffusive dynamics dominates 
the kinetics of phase transformations, the diffusion of vacancies, the motion of grain 
boundaries, and dislocation climb. In other words, elastic interactions mediated by 
wave modes have been incorporated, which travel on time scales that are orders of 
magnitude slower than the molecular vibrations yet considerably faster than the 
diffusive time scale. A similar EOM has been proposed for the VPFC model, whose 
free-energy functional forces the order parameter to be non-negative. The resulting 
approach dictates the number of atoms and describes the motion of each of them. 
Solution of the respective EOM might be viewed as essentially performing MD 
simulations on diffusive time scales [56]. A similar approach has been adapted to 
study the dynamics of monatomic and binary glasses, however, using MPFT-type 
free-energy functional and temperature-dependent noise terms [47] . 



3.1.1.3. The Euler-Lagrange equation. The ELE can be used to study equi- 
librium features including the mapping of the phase diagram [28] as well as the 
evaluation of the free energy of the liquid-solid interface [28] and of the nucleation 
barrier [35]. We note that noise may influence the phase diagram and other physical 
properties. Therefore, the results from the ELE and EOM are expected to converge 
for C ^ 0. We also call attention to the fact that so far as the equilibrium results 
(obtained by ELE) are concerned, the SH and IM-PFC models are equivalent. 
Once the free-energy functional is defined for the specific PFC model, its extremes 
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Figure 8. Multiplicity of localized hexagonal cluster solutions for the dynamic SH equation du/dt = 
— (1 + V|)^w — jJLU -\- uv? — with u = 1.6. The central panel displays a part of the bifurcation diagram 
for localized hexagonal patches (for details see reference [142]). The solid and dashed lines stand for the 
stable and unstable solutions, respectively. The vertical lines in grey correspond to the fold limits of planar 
(10) and (11) hexagon pulses. The red and green regions indicate where temporal self-completion does or 
does not occur, respectively. Panels 1-4 show colour plots of the hexagon patches at the inner and outer 
left folds. A different parametrisation is used here: fi = —e -\- SipQ. (Reproduced with permission from D. 
J. B. Lloyd, B. Sandstede, D. Avitable, and A. R. Champneys, Localized hexagon patterns of the planar 
Swift-Hohenberg equation, SIAM J. Appl. Dyn. Syst. 7 (2008), 1049-1100, no. 3, DOT: 10.1137/070707622 
© 2008 by The Society for Industrial and Applied Mathematics.) 



can be found by solving the respective ELE, which reads as 



(65) 



where '0o is the reduced particle number density of the unperturbed initial liquid, 
while a no- flux boundary condition is prescribed at the boundaries of the simulation 
window (n-Vf'0 = and (n-Vf)V?'0 = 0, where n is the normal vector of the 
boundary). For example, inserting the IM-PFC free energy and rearranging the 
terms, one arrives at 

(- e + (1 + V|f)i^P - V^o) - 4') - - V^o) = . (66) 

Equation (66) together with the boundary conditions represents a fourth-order 
boundary value problem (BVP). 

Multiplicity of solutions of the ELE. It is worth noting that in the case of the 
PFC/SH-type models, a multiphcity of solutions can usually be found for the same 
BVP, defined by the boundary conditions and the ELE. This feature of the sta- 
tionary solutions has been recently addressed in some detail in 2D for the SH [142] 
and in ID for the VPFC [143] models. (Figure 8 illustrates this phenomenon via 
showing the bifurcation diagram for compact hexagonal clusters in the SH model 
[142].) 

3.1.2. Binary PFC models 
3.1.2.1. The free energy. 



Binary Swift-Hohenberg form. The earliest binary extension of the PFC model 
has been proposed in the seminal paper of Elder et al. [12] obtained by adding an 
interaction term with coefficient A to the free energy of two mixed single-component 
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SH free energies. Working with particle densities ipi and V^2, this yields the free 
energy 




Here, the physical properties (bulk moduli, lattice constants, etc.) of the individual 
species are controlled by the parameters with subscripts 1 and 2, respectively, and 
by the average values of the particle densities ^1 and '02- It has been shown that 
the model can be used for studying structural phase transitions [12]. 

DFT-hased binary PFC model. The most extensively used binary generalisation 
of the IM-PFC model has been derived starting from a binary perturbative DFT, 
where the free energy is Taylor expanded relative to the liquid state denoted by 
and pb, respectively, up to second order in the density differences ApA = PA — PA 
and ApB = Pb — Pb (up to two-particle correlations) [26]: 

- \ Jdn jdr2 (ApA(ri)c^l(ri, r2)ApA(r2) ^^g^ 

+ ApB(ri)4B(ri,r2)ApB(r2) 
+ 2ApA(ri)4'^(ri,r2)ApB(r2)). 



It is assumed here that all two-point correlation functions are isotropic, i.e., 
qy(ri,r2) = q^- (|ri — r2|). Taylor expanding the direct correlation functions in 
Fourier space up to fourth order, one obtains c^^\\ri — r2|) = (4jO ~ 4j2^r2 + 

c[^\v^^)6{ri — r2) in real space [26]. The partial direct correlation functions c\^^ 
can be related to the measured or computed partial structure factors. 

Following Elder etal. [26], the reduced partial particle density differences are 
defined as ^pA = {pa-Pa)/po and = (pb-Pb)/po, where po = pa+Pb- They have 
introduced then the new variables ip = ipA + V^B and i(j = {ipB — '^a) + {pb — Pa) / Po^ 
obtaining fields whose amplitude expansion yields field amplitudes resembling order 
parameters associated with structural and composition changes in the conventional 
PF models. Expanding the free energy around lip = and ip = one obtains 



Pvef ^bT 



(69) 



This model has been used for studying a broad range of phase transitions, including 
the formation of solutal dendrites, eutectic structures [26, 28], and the Kirkendall 
effect [55]. 
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Binary PFC model for surface alloying. A somewhat different formulation of the 
binary model has been proposed to model compositional patterning in monolayer 
aggregates of binary metallic systems by Muralidharan and Haataja [48]: 

F = Jdr (|(/3(c) + + Vlf)p+^ + V{c)p 

+ fo (^(Vrc)2 - + ^ ((1 + c) ln(l + c) + (1 - c) ln(l - c)))) . 

(70) 

In this construction, the values c = ±1 of the concentration field stand for different 
atomic species, V{c) is the atom specific substrate-film interaction, kc incorporates 
the bulk lattice constant of different species, fo governs the relative strength of the 
elastic and chemical energies, while 0c is the critical temperature, the scaled ab- 
solute temperature, and wq tunes the chemical contribution to the interface energy 
between different species. Starting from this free-energy functional, Muralidharan 
and Haataja have shown, that their PFC model incorporates competing misfit dis- 
locations and alloying in a quantitative way, and then employed the model for 
investigating the misfit- and line tension dependence of the domain size [48] . 



3.1.2.2. The equations of motion. 



Dimensionless binary form in Swift- Hohenberg fashion. The same type of over- 
damped conservative dynamics has been assumed for the two types of atoms [12]: 



Here, Mi and Q are the mobility and the noise term applying to species i G {1, 2}. 



DFT-based binary PFC model. In this widely used formulation of the binary 
IM-PFC model, it is assumed that the same mobility M applies for the two species 
A and B (corresponding to substitutional diffusion) decoupling the dynamics of the 
fields ip and t/j. Assuming, furthermore, a constant effective mobility Mq = 2M/p^ 
and conserved dynamics, the equations of motions for the two fields have the form 
[26] 

^ = MeV?^ and ■i^M.Vl^. (72) 

In general, the coefficients Sg, and R in equation (69) depend on '0. A Taylor 
expansion of Bi^ip)^ ^5(^)5 and i?('^) in terms of ijj yields new coefficients S|, S|, 
and Ri with i = 0, 1, 2, . . . that are independent of ip. Retaining only the coefficients 
Bq, B2, Sq, i?o, and and inserting the free energy (69) into equations (72), one 
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obtains 
dtp 

~dt 



MeV? (^^'(^0 + B2lp^) + S^'^ + 

Bf 



dtp 

'dt 



+ ^(2{R^ + R,^)\l + {Ro + RitP^)vt)tl; (73) 

+ ^(2v2(v(i?o + Rii')^) + v^(v(i?o + Rii')'))) , 

MeV^ (^Bltptp"^ + 'y + tutp + utp^ - L^V'^tp 

+ 2B'otp({Ro + Ritp)RiVl + {Ro + RitpfRiVt)tp) ■ 



(74) 



3.1.2.3. The Euler- Lagrange equations. Since the derivation of the ELE is 
straightforward for all the binary PFC models, we illustrate it for the most fre- 
quently used version deduced from the binary Ramakrishnan-Yussouff-type clas- 
sical perturbative DPT. The extremum of the grand potential functional requires 
that its first functional derivatives are zero, i.e.. 



_ 

6ip 6ip 



and — - = — - 



(75) 



where tpo and t/^Q are the total and relative particle densities for the homogeneous 
initial state. Inserting equation (69) into equation (75), one obtains after rearrang- 
ing 



(76) 



(^B,tP + B.RtP^ (2V2 + i?^2y4^ + |£ (2V2(il2) + V^(i?4))j _ ^o) 

= -s{tp''-tpl)-v[tp''-tPl), 

L'WltP - ^{tptp^ - V^oVo') - 2B,R^tp{Vl + R'Vt)tP 

dip oy) (j'j^ 

where R = i?o + These equations are to be solved assuming no-flux boundary 
conditions at the border of the simulation box for both fields [n-Vr^ = 0, (n- 
Vr)V^V^ = 0, n-VrV^ = 0, and (n-Vr)V^V^ = 0]. 

3.1.3. PFC models for liquid crystals 

Another level of complexity is to consider interactions between particles, that 
are not any longer spherically symmetric. The simplest case are oriented particles, 
i. e., particles with a fixed orientation along a given direction where the interaction 
U2{yi —Y2) is not any longer radially symmetric [as assumed in equation (36)]. The 
derivation described in section 2.3 was straightforwardly extended towards this case 
leading to a microscopic justification of the APFC model (59) proposed in section 
3.1.1.1. Clearly, the resulting crystal lattices are anisotropic, which leads also to 
an anisotropic crystal growth. Much more complicated, however, is the case of 
orientahle particles, which can adjust their orientation and form liquid crystalline 
mesophases. 
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Therefore, we now address orientable anisotropic particles, that possess orien- 
tational degrees of freedom. In order to keep the effort manageable, the particles 
are assumed to be uniaxial, i.e., they are rotationally symmetric around an in- 
ternal axis. Consequently, their orientations are described by unit vectors with 
i G {1, . . . , N} along their internal symmetry axes. Most of the statistical quantities 
can suitably be generalized towards orientational degrees of freedom by including 
also the orientational configuration space on top of the translational configurational 
space. 



3.1.3.1. Statics. For orientable particles, the one-particle density is now defined 
via 

N 



p{r^u) = {Y,S{r-r,)6{u-u,)) (78) 



i=l 



and contains also the probability distribution of the orientations as expressed by 
the dependence on the orientational unit vector u. Its full configurational mean is 
given by the mean particle number density 

where the orientational integral denotes integration over the two-dimensional unit 
sphere S2. Now the external potential is ?7i(r, u) and couples also to the particle 
orientation. In general, also the pair-interaction potential f72(ri — r2,ui,U2) has 
an orientational dependence. The corresponding equilibrium Helmholtz free-energy 
functional ^[p(r, u)] can be split as usual in three contributions 

T[p{r, u)] = Jld[p(r, u)] + j;xc[p(r, u)] + j;xt[p(r, u)] (80) 

namely the ideal rotator-gas free-energy contribution 

J-id[p(r,u)] =kBT Jdr Jdu p{r,u){\n{A^p{r,u)) - l) , (81) 

the nontrivial excess free-energy functional J^exc[p(j^? u)], and the external free- 
energy contribution (see section 2.1) 

j;xt[/>(r,u)] = Jdr Jdu p(r,u)C/i(r,u) . (82) 

Analogously to equation (8), a functional Taylor expansion 

oo ^ 

Te.Mr,u)] = -Fi°l(Pref) + ^6^^^ -j;W[p(r,u)] (83) 

n=l 

with the nth-order contributions 

H^l[p{r,u)] = - jdv^..-jdv^ jdn^.--jdnn cH(r", u") nAp(r,, u^) , (84) 
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the n-particle direct correlation function c(^)(r",u"), and the abbreviations 

r- = (ri,...,rO, = (ui, . . . , (85) 

are used. Again, the constant zeroth-order contribution of the functional Tay- 
lor expansion (83) is irrelevant and the first-order contribution vanishes. In the 
Ramakrishnan-Yussouff approximation, the functional Taylor expansion is trun- 
cated at second order. 

In order to derive PFC-type models for liquid crystals, one can follow the same 
strategy as for spherical systems: first, the full density field p(r, u) is parametrized 
by small and slowly varying space-dependent multi-component order-parameter 
fields. An insertion into the density functional together with a perturbative and 
gradient expansion yields a local PFC-type free energy. Then several coupling terms 
of the order-parameter components arise, whose prefactors are given by moments 
of the generalized direct fluid correlation functions. Still, the gradient expansion is 
more tedious, since the orientational space has a more complicated topology. For 
stability reasons, one has to assume that the coefficients of the highest-order terms 
in the gradients and order-parameter fields in the PFC model are positive in the 
full free-energy functional. If this appears not to be the case for a certain system, 
it is necessary to take further terms of the respective order-parameter field up to 
the first stabilizing order into account. 



3.1.3.2. Two spatial dimensions. We now consider first the case of two spatial 
dimensions both for the translational and orientational degrees of freedom. Here, 
the topology of orientations is simpler than in three spatial dimensions. Obviously, 
all previous expressions can be changed towards two dimensions by replacing the 
volume V by an area A, by changing the unit sphere ^2 into the unit circle S'l, 
and by replacing A^ by A^ in equation (81). The orientational vector u{^) = 
(cos((/?), sin((/?)) can be parametrized by a single polar angle (p G [0,27r). 



Derivation of the PFC free- energy functional First we chose as order-parameter 
fields the reduced translational density^ the polarization^ and the nematic tensor 
field. The reduced translational density is defined via 

'0(r) = ydu (/9(r, u) - pref ) , (86) 

while the polarization is 

P(r) = (d\Jip{Y,\ji)\i (87) 

^Pref J 

and describes the local averaged dipolar orientation of the particles. Finally, the 
symmetric and traceless nematic tensor with the components 

Qij{Y) = /dup(r, u) (uiUj - l-6ij] (88) 



describes quadrupolar ordering of the particles. Equivalently, one could decompose 
the polarization P(r) = P(r)p(r) into its modulus P(r) and the local normalised 
dipolar orientation p(r) and use the two order parameters P{r) and p(r) instead 
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of P(r). Similarly, the nematic tensor can be expressed by [144, 145] 



Qi,{v) = S{v)(n,{v)n,{v)--5ij] (89) 



through the nematic order parameter S'(r), which measures the local degree 
of quadrupolar orientational order [146], and the nematic director n(r) = 
(ni(r), n2(r)). Notice that n(r) and p(r) do not necessarily coincide. 
One may expand the total density in terms of its orientational anisotropy as 

p(r,u) = p^Qi{l + ^l){Y) + Pi{Y)ui + UiQij{Y)uj) . (90) 

Inserting the parametrisation (90) into equation (81), performing a Taylor expan- 
sion of the integrand up to fourth order in the order-parameter fields, which guar- 
antees stability of the solutions, and carrying out the angular integration yields to 
the approximation 

^id[^, P^. Q^J] = i^id + TTPref k^T Jdv /id(r) (91) 

with the local scaled ideal rotator-gas free-energy density 



M = ^(8 - 2/>^ + 2P,Q,jPj - Ql) + ^(4 + 2/>^ + QD -\ + \ 



(92) 



where F^^ = 2n p^^i k^T A{\n{h? pj-Qf) — 1) is an irrelevant constant. 

If the functional Taylor expansion (83) for the excess free energy is truncated at 
fourth order, an insertion of the parametrisation (90) into the functional (84) and 
a gradient expansion yield [147] 

J-i^) Pi, Qij\ = -jdv fi-l (r) (93) 

with 

/i^) = AiV'" + A2{di^f + A^idl^f + Bi{d,^)Pi + B2Pi{djQij) + B:,{di^){djQij) 
+ C^Pf + C2Pi{dlPi) + C^idiPif + D^Q% + D2idjQijf , (94) 

fil = E^xl;^ + E2iPPl + E3i;Q^j + E^PQijPj + Fii?{diPi) + F2i;Pi{djQij) 
+ F3{M)QijPj + F^Pf{djPj) + F^{diPi)Qli + FePiQkiidjQkj) , 



(95) 
(96) 



Obviously, this creates a much more sophisticated series of coupling terms between 
gradients of the different order-parameter fields. They are all allowed by symmetry. 
There are altogether 28 coupling coefficients A^, B^, D^, Ei^ Fi^ G^, which can in 
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Table 2. Relevant special cases that are contained in the polar PFC model for two spatial dimensions, nmax is 
the order of the functional Taylor expansion (83). If nmax is not specified, it can be arbitrary (arb.). 
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Qij{r) 


2 


PFC model of Lowen [58] 


^(r) 




<5ii(r) 


4 


Full free-energy functional 



principle all be expressed as moments over the microscopic fluid direct correlation 
functions. These explicit expressions are summarized in appendix A.L 

The general result for the two-dimensional PFC free energy constituted by equa- 
tions (94)- (96) contains several special cases, known from the literature, which we 
now discuss in more detail. 

These special cases follow from the full free-energy functional either by choosing 
some of the order-parameter fields as zero or as a constant different from zero and 
by taking into account the contributions of the functional Taylor expansion (83) 
only up to a certain order nmax ^ {2,3,4}. Table 2 gives an overview about the 
most relevant special cases. 

The two most simple special cases are obtained for either a constant polarization 
Pi(r) or a constant nematic tensor Qijir) with an arbitrary choice for rimax, when 
all remaining order-parameter fields are assumed to be zero. These special cases are 
known to be Landau expansions in Pi{r) and Q^j(r), respectively. If the polarization 
Pi{r) is not constant, but space-dependent, and the other order-parameter fields 
vanish, then the functional can be called a gradient expansion in the polarization. 
Analogously, a gradient expansion in the nematic tensor Qij{r) is obtained, if 
it is space-dependent and all other order-parameter fields vanish. If additionally 
^max = 2 is chosen, this gradient expansion becomes the Landau-de Gennes free 
energy for inhomogeneous uniaxial nematics [145]. When only '^(r) is constant and 
the other order-parameter fields are space-dependent, we recover the case of an 
incompressible system. If all anisotropies are neglected, we recover with rimax = 2 
the original PFC model of Elder etal. [26]. The liquid crystalline PFC model of 
Lowen [58] is obtained, if in addition to the translational density also the nematic 
tensor is space-dependent and if it is parametrized according to equation (89). In 
the polar PFC model for two spatial dimensions, one can also consider the case of a 
space-dependent translational density '^(r), a space-dependent polarization Pi(r), 
and a vanishing nematic tensor, that corresponds to a ferroelectric phase without 
orientational order, but such a phase was never observed in experiments up to the 
present day. Therefore, this case is not included in Tab. 2. 

Equilibrium bulk phase diagram. While there are two independent parameters 
for the original PFC model (see equation (1) and figure 9), the number of coupling 
parameters explodes for the general liquid crystalline PFC model as proposed in 
equations (94)- (96) to 27 coupling coefficients. One of them can be incorporated 
into a length scale such that 26 coefficients remain. Therefore, a numerical explo- 
ration of the equilibrium phase diagram requires a much higher effort. Maybe not 
too surprising, such calculations are sparse and it was only until recently that the 
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phase diagram was calculated for the apolar PFC model [58], where P(r) = 0, 
which contains only 5 independent parameters [60]. We summarize and outline the 
basic findings of reference [60] in section 3.2.2. 



3.1.3.3. Three spatial dimensions. We now discuss the three-dimensional {d = 
3) PFC model for liquid crystals. In spherical coordinates, the three-dimensional 
orientation vector is 



u(0,0) = (sin(0)cos(0),sin(0)sin(0),cos(0)) (97) 

with the polar angle 6 G [0,7r] and the azimuthal angle (j) G [0, 27r). We consider 
here only the apolar case, where the polarization vanishes: P(r) = 0. There exists 
indeed a zoo of realizations of apolar particles both in the molecular and in the 
colloidal regime. For suitable interactions see references [148-156]. 

Following similar ideas as outlined in section 3.1.3.2 for two dimensions, we define 
as the corresponding order parameters the reduced translational density 

V^(r) = (p(r, u) - p^ef ) (98) 

and the 3 x 3-dimensional symmetric and traceless nematic tensor 

Qij{^) = ^^^^ Jdup{r, u) (^UiUj - ^6ij^ . (99) 

Again, the nematic tensor can be expressed by the nematic order-parameter field 
S{r) and the nematic director n(r) = (ni(r), n2(r), n3(r)), that is here the only 
unit vector that denotes a preferred direction in the liquid crystalline system. In 
the three-dimensional case, the decomposition of the nematic tensor is given by 
[144, 145] 

Q,,(r) = S{r) Qn,(r)n,(r) - ^S,j^ . (100) 

Notice that the nematic order-parameter field S{r) is the biggest eigenvalue of the 
nematic tensor Qij{r) and the nematic director n(r) is the corresponding eigenvec- 
tor. Accordingly, with the order-parameter fields V^(r) and Qij{r)^ the one-particle 
density is approximated by 

p(r,u) = pref(l + '^(r) +UiQij{r)uj) . (101) 

As before, the Helmholtz free-energy functional has to be approximated by a Taylor 
expansion around the homogeneous reference system and by a gradient expansion 
[59]. The ideal rotator-gas free-energy functional is approximated by 

^id[V^, Q^j] = i^id + vrpref ^rT Jdv /id(r) (102) 
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Table 3. Relevant special cases, that are contained in the apolar PFC model for three spatial dimensions. 
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with the local scaled ideal rotator-gas free-energy density 



15 315 y " V 15 y 3 3 

^ 4tr(Q^) _ 16tr(Q^) ^ 8tr(Q4) 



(103) 



15 315 315 



Here, tr( • ) denotes the trace operator and F^^ = AT{p^QiVk^T{\]i{h? p^^f) — 1) 
an irrelevant constant. For the excess free-energy functional, the Ramakrishnan- 
Yussouff approximation (11) is used together with equation (101) involving the 
direct correlation function c^'^\yi — r2,Ui,U2). Respecting all symmetries, one fi- 
nally obtains the following approximation for the excess free-energy functional 

4^[i^,Qij] = -Jdrf^^{r) (104) 

with the local scaled excess free-energy density 



+ KiidjQijf + K2Qij{dlQ. 



(105) 



The seven coupling coefficients Ai, ^2, A3, 5i, ^2, i^i, K2 can be expressed as 
generalized moments of the microscopic correlation function c^'^\ri — r2,ui,U2). 
The full expressions are summarized in appendix A. 2. We remark that Ki and K2 
correspond to the traditional Frank constants appearing in Frank's elastic energy 
[145, 157]. 

Special cases of the free-energy density are summarized in Tab. 3. As for the PFC 
model for two spatial dimensions, one obtains a Landau expansion in the nematic 
tensor Q^j(r), if the translational density ip{r) is zero and the nematic tensor is 
constant so that all gradients vanish. When only '^(r) is constant and Qij(r) is 
space-dependent, the Landau-de Gennes free energy for inhomogeneous uniaxial 
nematics [145] is recovered again. Clearly, the original PFC model of Elder et al. 
[26] for isotropic particles in three spatial dimensions can be obtained from the 
full free-energy functional by choosing Qijir) = 0. Finally, we note that the full 
functional was recently numerically evaluated for several situations by Yabunaka 
and Araki in reference [158]. 



3.1.3.4. Dynamics. Dynamical density functional theory (DDFT) can now be 
used to derive the dynamics of the PFC models for liquid crystals. DDFT is 
well justified for Brownian anisotropic particles (as, for example, colloidal rods 
or platelets). The basic derivation is similar to that performed for spherical parti- 
cles in section 2.3.2, but in practice it is much more tedious. The basic derivation 
is performed in three steps. At first, the order-parameter fields, that have been 
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chosen for the statics, are assumed to be time-dependent and the time-dependent 
one-particle density p(r, u, t) is approximated in terms of these time-dependent 
order-parameter fields. Secondly, the chain rule for functional differentiation is used 
to express the functional derivative 5J^ /5p of the Helmholtz free-energy functional 
T in terms of the functional derivatives of the free-energy functional with respect 
to the chosen order-parameter fields. Finally, the time-dependent parametrisation 
for the one-particle density and the time-dependent expression for the functional 
derivative 8T/5p are inserted into the DDFT equation and a set of in general 
coupled dynamic equations for the single order-parameter fields is obtained by an 
orthogonal projection of the DDFT equation with respect to the orientation u. 

We first consider the case of two spatial dimensions (c? = 2). The time-dependent 
noise- averaged one-particle number density is now 

p(r,u, t) = / ^5(r - r,(t))5(u - u,(t))\ (106) 
\ 1=1 ' 

and its order-parameter parametrisation is 

p(r,u,t) = p^^i{\ + ^(v,t) + Pi{v,t)ui + UiQi^{v,t)u^) . (107) 

The DDFT equation for orientational degrees of freedom in two spatial dimensions 
(without a hydrodynamic translational-rotational coupling) reads [99] 

. , _ /Dt(u) , . [p]) 

/)(r,U,t)=Vr- ^^/^(r.U,t)V: ^ '^^^^ 



. k^T ' ' ^p(r, u, i) , 

(108) 

+ ^5,^p(r,u,t)5,^-^^J 

with the translational short-time diffusion tensor 

Dt(u) = D\\Vi (g) u + Dx_(\ - u (g) u) . (109) 

Here, D|| and are the translational diffusion coefficients for translation parallel 
and perpendicular to the orientation u = (cos((/)), sin((/))), respectively, and the 
symbol 1 denotes the 2 x 2-dimensional unit matrix. The two terms on the right- 
hand-side of this DDFT equation correspond to pure translation and pure rotation, 
respectively. 

The constant-mobility approximation (CM A) is now 

p(r, U, t) - p,ef Vr • Vr^^^^ j + P-f ^ ^^^^^^ • (HO) 

Within the CMA, this yields for the dynamics of the order-parameter fields 

^ + dijf = 0, P, + $f = 0, Qij + ^f. = (111) 
with the currents and quasi-currents 



Sip SQij 



(112) 
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The coefficients in equations (112)-(114) are defined as 

L>ll + D_L £>|| + 3D_L D\\ - D± 

(115) 



Notice that Z^y ^ holds for all types of uniaxial particles, if the vector u for 
the orientation of the axis of symmetry is chosen properly. 

Finally, we remark that in three dimensions {d = 3) the dynamics is much more 
involved but can in principle be derived along similar lines from DDFT [103]. 
Recently, the CMA dynamics was considered in reference [158]. A main result of 
this reference is shown in section 5.3. 

3.1.4' Numerical methods 

3.1.4.1. The equation of motion. The commonly used explicit time-stepping 
with a finite-difference (FD) scheme has been routinely applied for solving the 
EOM of the PFC models [159]. However, owing to the high order differential op- 
erators appearing in the EOM of the PFC models (up to 12th order), explicit 
time-stepping suffers from severe constraints. Energy stable large time-step im- 
plicit FD methods have been developed for the PFC [160] and MPFC [161] equa- 
tions that lead to large sets of sparse algebraic equations. The resultant algebraic 
equations can be solved using nonlinear multigrid methods [162]. As elegant alter- 
native, pseudo-spectral methods can be used that combine unconditionally stable 
time marching while resulting in algebraic equations of the diagonal form [163]. 
Furthermore, pseudo-spectral methods offer exponential convergence with the spa- 
tial resolution as opposed to the polynomial convergence rate of the FD schemes. 
This means that with smaller spatial resolution one obtains results comparable to 
those obtained with high spatial resolution in the real space methods. Using opera- 
tor splitting techniques the domain of the spectral methods has been extended to a 
wide range of problems, including PDEs with variable coefficients [164]. A detailed 
study on the application of this method to the PFC models shows that the respec- 
tive schemes can be parallelised easily and efficiently [164], yielding up to 10^ times 
faster computations, when compared to FD schemes. Adaptive time stepping [163] 
may further accelerate the computations. Other numerical methods (e.g., stable 
semi-implicit finite element discretisation combined with adaptive time stepping 
[165]) have also been used for PFC models, however, further investigations are 
needed to assess their numerical accuracy. 



3.1.4.2. The Euler-Lagrange equation and other saddle point finding meth- 
ods. In recent works (see, e.g., references [28, 35]) the ELE has been solved us- 
ing a semi-spectral successive approximation scheme combined with the operator- 
splitting technique [166]. A different approach termed the fixed length simplified 
string method has been proposed to find the minimum energy path and the nucle- 
ation barrier in reference [32]. 
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3.1.5. Coarse- graining the PFC models 

Soon after the appearance of the IM-PFC model, attempts have been made to 
use this atomistic approach as a basis for deriving phase-field-type coarse-grained 
models using the amplitude equation method. This extension to the IM-PFC for- 
malism, when combined with the adaptive grid, has the potential to enable simula- 
tions of mesoscopic phenomena (/im mm) that are resolved down to the atomic 
scale, still incorporating all the respective physics. 

3.1.5.1. Amplitude equations based on renormalisation group theory. Gold- 

enfeld, Athreya, and Dantzig [131] have developed a computationally efficient ap- 
proach to polycrystalline solidification, based on the IM-PFC model. The nanoscale 
particle density distribution is reconstructed from its slowly varying amplitude and 
phase, obtained by solving the rotationally covariant equations of motion derived 
from renormalisation group theory. They have shown in two dimensions that the 
microscopic density distributions from their amplitude and phase equations show 
a very close match to the IM-PFC result. In later works Athreya and co-workers 
[19, 20] have combined this approach with adaptive mesh algorithms, leading to a 
substantial acceleration of the numerical code. (Possible ways of using the renor- 
malisation group methods have been discussed in references [20, 167].) 

3.1.5.2. Phenomenological amplitude equations. Approximate treatments 
based on the amplitude expansion of the free energy of the IM-PFC models (i.e., 
expressing in terms of the Fourier amplitude of the dominant density waves) have 
been developed. They have been used for various purposes, such as determining 
the anisotropy of liquid-solid interfacial free energies [41, 132], the Asaro-Tiller- 
Grinfeld (ATG) morphological instability of a stressed crystal surface, polycrys- 
talline growth from the melt, grain-boundary energies over a wide range of misori- 
entation, and grain-boundary motion coupled to shear deformation [54, 168]. Yeon 
et al. [169] have used an amplitude-equation approach to model the evolution of 
a two-phase system that has been validated by investigating the Gibbs-Thomson 
effect in two spatial dimensions. Elder, Huang, and Provatas [31] have proposed 
amplitude representations for the binary PFC model in the cases of triangular lat- 
tice (2D) and bcc and fee (3D) structures. The respective equations of motion have 
been related to those of the original PF theory of binary freezing and elasticity, 
providing explicit connection between the PF and PFC approaches. The abilities of 
the phenomenological amplitude models have been demonstrated for eutectic solid- 
ification, solute migration at grain boundaries, and for the formation of quantum 
dots on nanomembranes [31]. 

3.2. Phase diagrams the PFC models realize 

The different versions and extensions of the PFC model, reviewed in section 3.1, 
lead to different phase diagrams, which in turn depend on the dimensionality of 
the system. 

3.2.1. Phase diagram of single- component and binary systems 

The phase diagrams of the IM-PFC model [12] in two and three spatial dimen- 
sions are shown in figure 9. In two spatial dimensions, a single crystalline phase 
(the triangular phase) appears that coexists with the liquid and a striped phase 
[12]. Remarkably, phase diagrams of comparable structure have been predicted for 
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Figure 9. (a) Single-mode approximation to the phase diagrarn of the IM-PFC model in two spatial 
dimensions with e = (3 and the average reduced particle density ip. (Reproduced with permission from K. 
R. Elder, M. Katakowski, M. Haataja, and M. Grant, Modeling elasticity in crystal growth, Phys. Rev. 
Lett. 88 (2002), 245701, no. 24, DOI: 10. 1103/PhysRevLett. 88. 245701 © 2002 by the American Physical 
Society.) (b) Section of the 3D phase diagram of the IM-PFC model evaluated by the Euler-Lagrange 
method described in reference [28] with e = e and the average reduced particle density ipQ. Notice the 
stability domains of the bcc, hep, and fee phases. The homogeneous liquid is unstable right of the heavy 
grey line emerging from linear stability analysis [28]. (Reproduced from G. I. Toth, G. Tegze, T. Pusztai, 
G. Toth, and L. Granasy, Polymorphism, crystal nucleation and growth in the phase-field crystal model in 
2D and 3D, J. Phys.: Condens. Matter 22 (2010), 364101, no. 36, DOI: 10.1088/0953-8984/22/36/364101 
© 2010 by Institute of Physics Pubhshing.) 






(b) 



(c) 



Figure 10. (a) Temperature-particle density phase diagram from MD simulations with a DLVO-type 
potential. Circles represent the disordered phase, triangles the columnar phase, and squares the lamellar 
phase. Stars stand for points, where the free energies of two phases cross. Solid lines are a guide for the 
eyes. Snapshots of the (b) triangular rod and (c) lamellar phases. (Reproduced with permission from A. de 
Candia, E. Del Gado, A. Fierro, N. Sator, M. Tarzia, and A. Coniglio, Columnar and lamellar phases in 
attractive colloidal systems, Phys. Rev. E 74 (2006), 010403R, no. 1, DOI: 10.1103/PhysRevE.74.010403 
© 2006 by the American Physical Society.) 



weakly charged colloids with competing interactions [170]. In three spatial dimen- 
sions, as implied earlier by EOM studies [42], and confirmed by full thermodynamic 
optimization [171] and an equivalent method based on solving the ELE [28], stabil- 
ity domains exist for the homogeneous fluid, bcc, fee, and hep structures, besides 
the 3D version of the respective 2D structures: the rod and the lamellar structures. 
Interestingly, the rod and lamellar structures, and a phase diagram resembling to 
the IM-PFC phase diagram appear in MD simulations for a Derjaguin-Landau- 
Verwey-Overbeek (DLVO) type potential [172] (see figure 10), characteristic to 
charged colloidal systems. The IM-PFC model prefers the formation of the bcc 
phase near the critical point. 

In contrast, the 2M-PFC model of Wu et al. [132] designed to realize fee crystal- 
lization, suppresses the bcc phase [see figure 11(a)]. The 2M-PFC extension incor- 
porates the IM-PFC model as a limiting case. Interpolation between IM-PFC and 
the full fee limit {Ri = 0), in terms of the parameter Ri leads to the appearance of 
a bcc stability domain in the vicinity of the critical point [see figure 11(b)]. {Ri is 
the ratio of the Fourier amplitudes for the density waves having the second and first 
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(a) (b) 

Figure 11. (a) Single-mode approximations to the phase diagram of the 2M-PFC model in 3D for Ri = 0. 
(b) The same for Ri = 0.05. Notice the small bcc stability domain near the critical point. (Reproduced with 
permission from K.-A. Wu, A. Adland, and A. Karma, Phase- field- crystal model for fee ordering^ Phys. 
Rev. E 81 (2010), 061601, no. 6, DOI: 10.1103/PhysRevE.81. 061601 © 2010 by the American Physical 
Society.) 
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Figure 12. Free-energy density for the liquid (solid) and bcc (dashed) Fe as a function of reduced particle 
density in the EOF-PFC model. Crosses denote the equilibrium points obtained by the common tangent 
method. (Reproduced with permission from A. Jaatinen, C. V. Achim, K. R. Elder, and T. Ala-Nissila, 
Thermodynamies of bee metals in phase- field- erystal models^ Phys. Rev. E 80 (2009), 031602, no. 3, DOI: 
10. 1103/PhysRevE.80. 031602 © 2009 by the American Physical Society.) 

neighbour RLVs as wave vector.) Whether here the appearance of the bcc phase 
is accompanied with the appearance of a hep stabihty domain, as in the IM-PFC 
hmit, requires further investigation. 

The apphcabihty of the EOF-PFC model has been demonstrated for Fe [133]. 
The free energy vs. particle density curves for the solid and liquid phases, which 
were used to determine the equilibrium conditions at the melting point, are shown 
in figure 12. No phase diagram has been published for this model. It appears that 
much like the original, the IM-PFC model, it prefers bcc freezing. 

An attempt to control the preferred crystal structure relies on manipulating the 
two-particle direct correlation function so that its peaks prefer the desired struc- 
tural correlations [135]. The type of phase diagram accessible for this method and 
the respective free-energy curves for coexistence between the bcc and fee structures 
are displayed in figure 13(a). Starting from the observation that the nonlinearities 
can stabilize the square lattice [173] so that it coexists with the triangular phase 
[174], Wu, Plapp, and Voorhees [136] have proposed a method based on nonlinear 
resonance for constructing PFC models that prefer the desired crystal structure, 
as they indeed demonstrated for the square lattice in two spatial dimensions: the 
relative free energies of possible competing structures are compared in figure 13(b). 

Elder etal. [26] have obtained eutectic phase diagrams within the binary gen- 
eralisation of the IM-PFC model in two spatial dimensions both numerically and 



July 3, 2012 0:16 Advances in Physics aip 

Advances in Physics 41 




n \|/ 

(a) (b) 

Figure 13. (a) Phase diagram in the effective temperature cr-mean density n-plane for a system whose di- 
rect correlation function was manipulated so that peritectic coexistence between the fee and bcc structures 
is realized. (Reproduced with permission from M. Greenwood, J. Rottler, and N. Provatas, Phase-field- 
crystal methodology for modeling of structural transformations, Phys. Rev. E 83 (2011), 031601, no. 3, 
DOT: 10. 1103/PhysRevE.83. 031601 © 2011 by the American Physical Society.) (b) Free-energy density 
/ vs. average reduced particle density ip when controlling the crystal structure via nonlinear resonance. 
/l denotes the free-energy density of the liquid. The analytical solutions are plotted as lines, while the 
numerical simulation results for rolls, squares, hex-1, and hex-2 are plotted as squares, circles, diamonds, 
and triangles, respectively. (Reproduced with permission from K.-A. Wu, M. Plapp, and P. W. Voorhees, 
Controlling crystal symmetries in phase-field crystal models, J. Phys.: Condens. Matter 22 (2010), 364102, 
no. 36, DOT: 10.1088/0953-8984/22/36/364102 © 2010 by Institute of Physics Publishing.) 




(a) (b) (c) 

Figure 14. (a), (b) Reduced temperature ABq vs. reduced particle density difference ip phase diagrams for 
the two-dimensional triangular system in the amplitude equation formalism. The filled regions correspond 
to regions of phase coexistence. (Reproduced with permission from K. R. Elder, Z.-F. Huang, and N. 
Provatas, Amplitude expansion of the binary phase-field-crystal model, Phys. Rev. E 81 (2010), 011602, 
no. 1, DOT: 10. 1103/PhysRevE.81. 011602 © 2010 by the American Physical Society.) (c) Colour map for 
the thermodynamic driving force of eutectic solidification (Acj is the grand potential density difference 
relative to the liquid) for the binary IM-PFC model in 3D [28], as a function of the chemical composition 
ipo and density ipo. Notice that solidification is expected in the region, where Acu < 0. (Reproduced from 
G. I. Toth, G. Tegze, T. Pusztai, G. Toth, and L. Granasy, Polymorphism, crystal nucleation and growth 
in the phase-field crystal model in 2D and 3D, J. Phys.: Condens. Matter 22 (2010), 364101, no. 36, DOT: 
10.1088/0953-8984/22/36/364101 © 2010 by Institute of Physics Publishing.) 

by using the single-mode approximation. Comparable phase diagrams [see figures 
14(a) and 14(b)] have been reported for triangular phases in two spatial dimen- 
sions, and for the bcc and fee phases in 3D, by the phenomenological amplitude 
equation method [31]. The 3D extension of the binary IM-PFC model has been 
investigated by Toth etal. [28]. The map of the thermodynamic driving force for 
solidification as a function of composition and density of the initial liquid is shown 
in figure 14(c). 3D eutectic solidification to bcc phases of different lattice constants 
has indeed been observed in the domain of the largest driving force. 

An essential question is how to extend the PFC framework to accommodate 
the physical properties of real alloys of different crystalline structures. A two- 
phase binary extension of the EOF-PFC model based on the IM-PFC and 2M-PFC 
concepts might be worth exploring. 

The VPFC extension of the IM-PFC model leads to a restructured phase diagram 
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Figure 15. Phase diagram for the VPFC model in two spatial dimensions [see equation (58)] (a) for 
ko = 1 [143]. The red solid lines denote coexistence curves, the green dash-dotted lines envelope the 
region, where localized and hexagonally ordered density peaks (bumps) coexist, whereas the blue dashed 
line indicates the linear stability limit of the spatially uniform phase. Snapshots of simulations for (b) 
stripes and (c) hexagonally ordered holes are also shown. These simulations have been performed for: 
ko = 1, p = 0.9, h = 1500, M = 1 (the mobility in the EOM) and (b) ipo = 0.4 and (c) ipo = 0.53. 
(Reproduced from M. J. Robbins, A. J. Archer, U. Thiele, and E. Knobloch, Modelling fluids and crystals 
using a two- component modified phase field crystal model, Phys. Rev. E 85 (2012), 061408, no. 6, DOI: 
10. 1103/PhysRevE.85. 061408 © 2012 by the American Physical Society.) 

(see figure 15) [143], whose central region contains stability domains for the 2D 
hexagonal crystal and localized density peaks ( "bumps" ) that represent individual 
particles. 

3.2.2. Phase diagram of two-dimensional liquid crystals 

For liquid crystalline systems, there are much more candidates for possible bulk 
phases. Therefore the topology of the bulk phase diagram is getting more complex. 
Recently, bulk phase diagrams were computed by using the two-dimensional apolar 
PFC model of Lowen [58]. After an appropriate scaling in energy and length, this 
model reads as 

^[^, Qi,] ^Jdr - + _ + - 1)^ + 

+ A^xl? + A^i>{dl + dldf)i^ (116) 
+ B^{diil,){djQij) + D^Ql + D2{djQijf) . 

Some trends of the phase diagram can directly be read of equation (116). Since the 
parameter Di controls the contribution of the nematic tensor Qij (r) and therefore 
also of the nematic order parameter S'(r), the nematic phase can be expected to be 
stable for large negative values of Di. In the opposite case, if Di is large enough 
and positive, the term DiQfj(r) + (5|^(r)(5^^(r)/64 dominates the free energy and 
only phases with Qijiv) oc S{r) = can be stable. Crystalline phases with a non- 
vanishing nematic order can therefore only appear in a region around Di = 0. 
From previous work it is known that the difference Ai — ^2/4 has a big influence 
on the translational density field ^(r) [26, 60]. If the parameter Ai is large and 
positive, variations of the translational density field enlarge the free energy. Sim- 
ilarly, gradients of the translational density field enlarge the free energy for large 
and negative values of A2. Therefore, phases without any density modulations, i. e.. 
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Figure 16. Stable liquid crystalline phases. The contour plots show the order-parameter fields '^(r) and 
»S'(r) in the a;i-a;2-plane for the isotropic and nematic phases, the stripe phase and columnar/smectic A 
phase, two plastic triangular crystals with different orientational ordering, and a plastic honeycomb crystal 
as well as a plastic square crystal. The black lines in the plots of the second and fourth column represent the 
director field n(r). In the plots with S(r) = 0, the director field is not shown because it is not defined. The 
parameters are A2 = 14, D2 = 8, and Bs = for the stripe phase and the plastic triangular crystal 1 and 
A2 = 14, D2 = 0.8, and Bs = —4 for all other phases. (Reproduced from C. V. Achim, R. Wittkowski, and 
H. Lowen, Stability of liquid crystalline phases in the phase- field- crystal models Phys. Rev. E 83 (2011), 
061712, no. 6, DOI: 10. 1103/PhysRevE.83. 061712 © 2011 by the American Physical Society.) 



the isotropic and the nematic phase, are preferred for positive values of the differ- 
ence Ai — A2/4, while all other phases with a periodic translational density field 
are preferred for negative values of this difference. Furthermore, there is a sym- 
metry concerning the reversal of the sign of the parameter in the free-energy 
functional. From equation (116) follows directly, that the free-energy functional is 
invariant under a simultaneous change of the signs of the parameter and the 
nematic order-parameter field ^(r). Due to this symmetry, S3 ^ can be assumed 
without loss of generality. 

Depending on the coupling parameters, a wealth of different stable liquid crys- 
talline phases was found (see figure 16). They include an isotropic phase, which 
has no translational and no orientational ordering ('0(r) = and S{r) = 0), a 
nematic phase with pure orientational ordering ('^(r) = and S(r) > 0), a colum- 
nar phase (or equivalently a smectic A phase), where the translational density and 
the orientation show a one-dimensional undulation, and various plastic crystals. In 
the columnar/smectic A phase the system has positional ordering in one direction, 
while it is isotropic perpendicular to this direction. The nematic order-parameter 
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Figure 17. Topological defects in three different plastic liquid crystals in the a::i-a::2-plane (schematic). The 
defects coincide with the maxima (red discs) and minima (cyan discs) of the translational density field '^(r). 
The symbols in the plots represent the following defects: (a) vortices with the topological winding number 
m = 1, (b) disclinations with m = — 1/2, (c) sources/sinks with m = 1, and (d) hyperbolic points with m = 
— 1. (Reproduced from C. V. Achim, R. Wittkowski, and H. Lowen, Stability of liquid crystalline phases in 
the phase- field- crystal model, Phys. Rev. E 83 (2011), 061712, no. 6, DOI: 10. 1103/PhysRevE.83. 061712 
© 2011 by the American Physical Society.) 



field S{t) for this piiase has a similar profile to the reduced translational density 
field '^(r) with maxima of these two fields at the same positions. Near the maxima 
of the translational density ^(r), the director field n(r) is preferentially parallel 
to the gradient 5^^(r), while it is perpendicular to 5i'0(r) around the minima of 
'^(r). A similar flipping of the orient ational field from perpendicular to parallel 
to the stripe direction was identified as transverse intralayer order in the three- 
dimensional smectic A phase of hard spherocylinders [109]. 

Plastic crystals are two-dimensional modulations of the translational density 
and have no global (averaged) orientational direction. Interestingly, depending on 
the coupling-parameter combinations, three different types of plastic crystals are 
stable including two-dimensional triangular, square, and honeycomb lattices for the 
translational ordering. The orientational field of plastic crystals is schematically 
shown in figure 17. It exhibits an interesting defect structure in the orientation. 
There are defects at the density peaks and in the interstitial places of the lattice. 
The defect structure has not yet been explored so far outside the PFC-world. 
It can be confirmed in microscopic computer simulations [175] or in real-space 
experiments [176-178]. 

In a two-dimensional slice of the coupling parameter space, the phase diagram 
is shown in figure 18. Apart from the isotropic phase with '^(r) = and S{r) = 0, 
which appears for Ai > A2/4 and Di > 0, for negative and large Di a nematic 
phase is stable. A rich topology occurs around i?i = including columnar phases 
and three different plastic crystals. However, one should bear in mind that the 
phase transitions shown in figure 18 were assumed to be isochoric, i. e., without any 
density jump. While this is in general a good approximation for liquid crystalline 
phase transitions, it is not true in general. Independent of the particular value 
of the parameter B^, the phase transition between the isotropic and the nematic 
phase is continuous, while all other phase transitions are discontinuous. 

In conclusion, rich phase diagrams with novel liquid crystalline phases were found 
in the apolar PFC model. This gives confidence that in a further step interfaces 
between coexisting phases and the dynamics can be explored on the basis of the 
PFC approach. It would open the way of PFC models to enter into the rich world 
of liquid crystals. 
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Figure 18. Phase diagrams for the parameters A2 = 14 and D2 = 0.8. The relevant liquid crystalline 
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resolution of the parameter space. (Reproduced from C. V. Achim, R. Wittkowski, and H. Lowen, Stability 
of liquid crystalline phases in the phase- field- crystal model, Phys. Rev. E 83 (2011), 061712, no. 6, DOI: 
10. 1103/PhysRevE.83. 061712 © 2011 by the American Physical Society.) 



3.3. Anisotropies in the PFC models 

One of the most attractive features of the PFC-type models is that the anisotropies 
for various physical properties follow directly form the crystal structure. The 
anisotropy of the interfacial free energy has been addressed theoretically and nu- 
merically by several authors [33, 40-44], whereas the growth anisotropy has been 
evaluated numerically in 3D [42]. 

3.3.1. Free energy of the liquid- solid interface 

3.3.1.1. Numerical results. Backofen and Voigt [43] have determined the 
anisotropy of the crystal- liquid interfacial free energy for small 2D clusters from 
simulations performed using the IM-PFC model. They have observed a strong 
dependence of the interfacial free energy on the distance from the critical point. 
The results have been fitted with the formula of Stashevich etal. [179] from low 
temperature expansion. Remarkably, the anisotropy shows strong size-dependence, 
when reducing the cluster size to a few particles. Comparable results have been 
obtained by Granasy etal. [33, 38] for the equilibrium shapes in a broader reduced 
temperature range (see figure 19). Toth etal. [28] evaluated the free energy and 
thickness of flat interfaces as a function of e via solving the ELE for the equi- 
librium (101) triangular crystal-liquid interface, and have shown that mean-field 
critical exponents apply. A similar approach has been applied in 3D for a large 
number of orientations of the flat bcc-liquid interface forming at e = 0.3748. The 
orientation dependence of the interfacial free energy and the respective Wulff plot 
are presented in figure 20 [44]. Apparently, the rhombic-dodecahedral equilibrium 
shape [42] obtained from simulations performed using the EOM for the same e has 
been a growth form. It became also evident that the usual cubic harmonic expan- 
sion is not sufficient for reproducing reasonably the anisotropy even if terms up to 
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Figure 19. Equilibrium shape for the IM-PFC model as determined by solving the EOM vs. reduced 
temperature e in the absence of noise [33]. (a) - (h): e = 0.05,0.10,0.15,0.20,0.25,0.30,0.325, and 0.35. 
Notice that the interface thickness decreases while the anisotropy increases with an increasing distance 
from the critical point. The computations have been performed on a 1024 x 1024 rectangular grid (the 
upper right quarter of the simulations is shown), whereas the crystalline fraction was ^ 0.3. Equilibra- 
tion has been performed for a period of 106 dimensionless time steps. Reduced particle density maps are 
shown. (Reproduced from L. Granasy, G. Tegze, G. I. Toth, and T. Pusztai, Phase-field crystal mod- 
elling of crystal nucleation, heteroepitaxy and patterning, Philos. Mag. 91 (2011), 123-149, no. 1, DOI: 
10.1080/14786435.2010.487476 © 2011 by Taylor &; Francis.) 




(a) 



(b) 



(c) 



Figure 20. Anisotropy of the bcc-liquid interfacial free energy at e = 0.3748 evaluated numerically for 
the single-component IM-PFC model and the respective equilibrium shapes [44]. (a) Gamma plot. Red 
dots represent the directions for which the interfacial free energy has been evaluated. The surface is a 
seventh-order cubic harmonics fit. Heavy dots are on the viewer's side of the surface. Light dots are on the 
other side of the surface, (b) Wulff shape evaluated from the red dots in the gamma plot, (c) The Wulff 
shape evaluated from a seventh-order cubic harmonics fit to the gamma plot. 

seventh order are considered. 



3.3.1.2. Analytical results. Wu and Karma [40] have used multi-scale analysis 
close to the critical point (0 < e <C 1) to evaluate the anisotropy of the interfacial 
fee energy. They have approximated the EOM of the IM-PFC model by a set of 
coupled equations describing the time evolution of the amplitudes of the dominant 
density waves. The analysis of the stationary solution led to an anisotropy that is 
independent of the reduced temperature e. This finding accords with those of Ma- 
janiemi and Provatas [41], who have used a simple coarse-graining technique, the 
local volume averaging method, for deriving amplitude equations for liquid-solid 
interfaces broad relative to the periodicity of the crystalline phase. In both studies, 
the temperature-independent anisotropy is a direct consequence of the approxima- 
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tions that lead to weakly fourth-order amplitude theories of the Ginzburg-Landau 
type, in which all material parameters can be scaled out from the free-energy 
functional [40, 41]. Accordingly, the anisotropy of the interfacial free energy de- 
pends only on the crystal structure. However, this independence of the anisotropy 
from temperature is unphysical, as it must vanish when the correlation length (the 
width of the liquid-solid interface) diverges in the critical point. In contrast to these 
phenomenological coarse-graining techniques, the renormalisation-group-based ap- 
proaches by Athreya, Goldenfeld, and Dantzig [131] lead to higher-order amplitude 
equations, from which the temperature cannot be scaled out; i.e., a temperature- 
dependent anisotropy is expected as found in the virtually exact numerical studies. 

3.3.2. Growth anisotropy 

Tegze etal. [42] have investigated the growth rate anisotropy in the IM-PFC 
model at e = 0.3748 for freezing to the bcc, hep, and fee structures. Diffusion- 
controlled layerwise crystal growth has been observed, a mechanism that is con- 
sistent with two-dimensional nucleation. The predicted growth anisotropics were 
found to decrease with increasing thermodynamic driving force (or velocity) con- 
sistently with kinetic roughening expected on the basis of IM-PFC simulations 
performed in two spatial dimensions [38] . 



3.4. Glass formation 

One of the most intriguing phenomena that may happen during solidification is 
glass formation, a process by which the undercooled liquid is transformed into an 
amorphous solid. An early IM-PFC study for single-component system by Berry 
etal. [46] relying on conservative overdamped dynamics indicated a first-order 
phase transition for this phase change, whereas the amorphous structure resembled 
closely to the glass structure obtained by embedded-atom-potential (EAP) molec- 
ular dynamics simulations for glassy Fe or Ni [180, 181]. These findings have been 
confirmed for the IM-PFC model by Toth etal. [140] and for the EOF-PFC model 
fitted to Fe by Toth et al. [28] and Granasy et al. [33]. In a more recent study. Berry, 
and Grant [47] have addressed glass formation in the framework of the monatomic 
and binary versions of the VPFC model with equation(s) of motion like equation 
(64) considering inertia and damping. They have shown that important aspects of 
glass formation can be reproduced over multiple time scales, including the agree- 
ment with mode-coupling theory (MCT) for underdamped liquids at low under- 
coolings and a rapidly growing dynamic correlation length that can be associated 
with a fragile behaviour. It appears that in the original PFC model glass formation 
takes place via a first-order phase transition, while the monatomic VPFC model 
behaves like poor glass formers, whereas the binary VPFC model displays features 
consistent with good glass formers. 



3.5. Phase- field- crystal modelling of foams 

Working at extreme distances from the critical point, Guttenberg, Goldenfeld, and 
Dantzig [61] have shown that the IM-PFC model can be used to describe the 
formation of foams. Under such conditions, the free energy of the periodic phases 
exceeds that for a mixture of two immiscible liquids - a situation that leads to 
the appearance of a foam-like structure coarsening with time. Starting from this 
observation, the authors present a simple PFC-type continuum scalar theory of wet 
and dry foams (see figure 21). 
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Figure 21. Foam structures predicted by the MPFC model, (a) Coexisting atoms and large cell foam, (b) 
Dry foam with no residual atoms, (c) Wet foam with circular bubbles. (Reproduced with permission from 
N. Guttenberg, N. Goldenfeld, and J. Dantzig, Emergence of foams from the breakdown of the phase field 
crystal model, Phys. Rev. E 81 (2010), 065301R, no. 6, DOI: 10.1103/PhysRevE.81.065301 © 2010 by the 
American Physical Society.) 

3.6. Coupling to hydrodynamics 

An interesting approach has been brought forward recently by Ramos etal. in 
reference [182]. It couples the PFC dynamics of the form of equation (1) to a 
Stokes-type dynamics as well as to an external pinning potential. The general set 
of equations the authors arrive at is given by 

dp ^ 

~dt ^ ~^ ~5p~ ~ ' 

where the specific configurational energy contribution i^int needs to be specified 
depending on the system. Assuming a periodic system, it is specified in reference 
[182] as 

i/int = ydr(^(-6+(l + V^)')V^ + ^+V^C/p). (118) 

In equations (117), p{r) and g(r) refer to particle and momentum density, re- 
spectively. Moreover, fi denotes the components of an external force vector and 
z/^(r,t) denotes noise components. Furthermore, 77 is a dissipative coefficient and 
Up accounts for the pinning potential. 

Even if at first glance the coupled set of equations (117) looks hke a first step 
towards coupling of the PFC model based on equation (1) to a full Navier-Stokes 
(NS) type equation and thus to hydrodynamic motion, a respective extension of 
equations (117) is not straightforward due to the simple relation between particle 
density p(r) and momentum density g(r) assumed. Thus instead of providing a 
step towards an extension of PFC models towards hydrodynamics, equations (117) 
rather provide a general framework for the inclusion of inertia effects in PFC mod- 
els. The latter issue has also been addressed in references [128, 183, 184], where it 
was shown that the consideration of inertia terms allows to include fast degrees of 
freedom into the PFC approach. Furthermore, the authors of reference [128] have 
shown that inclusion of such fast degrees of freedom yields a two-stage relaxation 
process of the system. Chen et al. proved the thermodynamic consistency of such 
models with fast degrees of freedom [185]. Thereby, these authors validated from 
the viewpoint of thermodynamic theory that such a two-stage relaxation process 
can truly be regarded as an important qualitative feature of nonequilibrium pattern 
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formation in periodic systems. To do so, they developed a unified thermodynamic 
framework applying to both conventional PF models and PFC models with slow 
and fast degrees of freedom. Based on this framework it is now possible to validate 
also PFC models based on their thermodynamic consistency, as done previously 
in numerous cases for PF models (see reference [186] for examples). Experimental 
evidence of inertia contributions to the dynamics of a non- linear system in the form 
described by the PFC model with fast degrees of freedom can indeed be found for 
example in rapid solidification [128]. 

The coupling of PFC models to hydrodynamics, though, is still an open issue. It 
makes sense to distinguish two cases: 

Case (a): Liquid-solid systems as, e.g., solidifying alloys, where hydrodynamic 
transport takes place in only the liquid phase. This implies that such systems can 
be treated in a multi-scale approach, where the PFC equation for the atomic scale 
dynamics in the interior of the solid is coupled to the long-range transport processes 
in the exterior, which - to be numerically efficient - couple only to the amplitude 
equation of the PFC equation. The reason is that in the case of PFC models it is the 
amplitude equation, which distinguishes between solid and liquid. Such amplitude 
equation approaches have been derived and solved numerically as computationally 
efficient representatives of PFC models, e.g., for diffusion limited poly crystalline 
grain growth [20] based on renormalisation group techniques (see section 3.1.5.1 
for details). Since renormalisation for PFC models of equation (1) type does not 
directly yield NS-type equations and thus a coupling to hydrodynamics exploiting 
the above multi-scale idea, this is still an open issue. 

Case (b) differs from the above picture in that it applies to systems, for which in- 
teratomic respectively intermolecular fluidic motion needs to be taken into account 
as, e. g., in polymers. Then one would desire a PFC model, where the PFC equation 
itself is directly coupled to a NS-type equation. A first step in that direction was 
done by Praetorius etal. [187]. They developed a PFC model with an advective 
term to simulate particles in a flowing solvent. To do so they followed the same 
route as Rauscher etal. [188] and Penna etal. [189] to derive a DDFT model for 
such systems, and approximated the latter further based on the Ramakrishnan- 
Yussouff approximation [85]. Strictly speaking the resulting PFC model applies 
only to potential flows. Furthermore, the coupling of the PFC variable ip into the 
dynamic equations for the hydrodynamic field is constructed simply based on nu- 
merical arguments. More rigorous generalisations of these studies might require 
additional physically motivated coupling terms, as demonstrated for hydrogels in 
terms of the derivation of the corresponding PF model in reference [190] (see section 
5.3 for details). 

Just as in the PF case, one can expect efficient further advance of PFC models 
in the future. Progress will result from physical modelling concepts based upon 
the classical routes to patterned non-linear systems [1], analytical mathematical 
expansion techniques as the renormalisation group approach (see section 3.1.5.1 
for details), and advanced numerical approaches [187, 191] in close comparison to 
experimental techniques. 

4. Phase-field-crystal models applied to nucleation and pattern formation in 
metals 

Crystal nucleation can be handled in two different ways within the framework of 
PFC models [28, 33, 35]: (i) via finding the properties of the critical fluctuations 
(nuclei) by solving the ELE under appropriate boundary conditions (zero field 
gradients at the centre and unperturbed liquid phase in the far field), whose so- 
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Figure 22. Particle density profiles for the critical fluctuations forming at different supersaturations (su- 
persaturation decreases from a to f), as obtained by the adaptation of the string method for determining 
the extrema of the free-energy functional in the IM-PFC model in two spatial dimensions. (Reproduced 
with permission from R. Backofen and A. Voigt, A phase- field- crystal approach to critical nuclei^ J. Phys.: 
Condens. Matter 22 (2010), 364104, no. 36, DOI: 10.1088/0953-8984/22/36/364104 © 2010 by Institute 
of Physics Publishing.) 

lution represents an extremum of the free-energy functional; (ii) by adding noise 
to the EOM. These approaches have their hmitations. (i) is expected to work for 
small undercoolings, where the individual heterophase fluctuations do not inter- 
act. Furthermore, it is not immediately straightforward how one should address 
possible non-crystalline nucleation precursors. In turn, in the case of (ii) it is not 
clear conceptually, which fraction of the thermal fluctuations is already integrated 
into the free energy and which wavelengths should yet be added as noise to the 
EOM [95, 192, 193] - a question inherently related to the proper choice of the 
high frequency cutoff one needs to make to avoid an ultraviolet catastrophe in 3D 
[194, 195]. Furthermore, the addition of noise to the EOM changes the free energy, 
the phase diagram, and the interfacial properties. While, in principle, correction of 
these is possible via parameter renormalisation [196, 197], further study is needed 
in the case of PFC models. On the other hand, the original free-energy functional 
used in (i) seems to miss the effect of longer wavelength fluctuations, which could 
move the system out of a metastable state. Considering them, (i) and (ii) may be 
regarded as approaches that provide complementary, probably qualitative informa- 
tion of the crystallizing system, which converge when the amplitude of the noise 
tends to zero. 

We note, however, that the results obtained via route (i) are more general than 
those obtained from (ii), as they follow directly from the free-energy functional, 
being therefore independent of the type of dynamics the EOM defines. Accordingly, 
results from (i) are valid even in cases, where the colloid- type diffusive dynamics is 
not applicable (e.g., metals). In this section, we review results obtained following 
route (i) for metals, or when pattern formation is governed by either chemical or 
surface diffusion. In turn, results obtained with the single-component PFC models 
with overdamped diffusive dynamics will be reviewed in section 4.3. 

4.1. Properties of nuclei from extremum principles 

4.1.1. Homogeneous nucleation 

An adaptation of the string method to find the saddle point of the free-energy 
functional has been used by Backofen and Voigt [32] for determining the proper- 
ties of the critical fluctuations in the IM-PFC model in two spatial dimensions. 
The respective density profiles are shown in figure 22. It is evident that at large 
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Figure 23. Nucleation barrier vs. size relationship obtained by solving the ELE for faceted nuclei (e = 
0.5) for the IM-PFC model [28] in two spatial dimensions: (a) homogeneous crystal nuclei for ipQ = 
—0.5134 + 0.0134/2^, where n = 0,1, 2,..., 7, respectively, (b) Heterogeneous crystal nuclei for = 
—0.5139 + 0.002/2^, where n = 0, 1, 2, . . . , 7, respectively. The lattice constant of the substrate is equal 
to the interparticle distance in the triangular crystal. Notice the substantial reduction of the nucleation 
barrier, a monolayer adsorption layer, and the contact angle of 60° determined by the crystal structure. 
(Reproduced from G. I. Toth, G. Tegze, T. Pusztai, G. Toth, and L. Granasy, Polymorphism, crystal 
nucleation and growth in the phase-field crystal model in 2D and 3D, J. Phys.: Gondens. Matter 22 (2010), 
364101, no. 36, DDI: 10.1088/0953-8984/22/36/364101 © 2010 by Institute of Physics Publishing.) 

supersaturations there are no bulk crystal properties at the centre of the smallest 
nuclei. They have also reported that small nuclei are faceted even though the large 
crystals are not. 

Toth et al. [28] have solved in two spatial dimensions the ELE of the IM-PFC 
model for the appropriate boundary conditions (homogeneous supersaturated liq- 
uid in the far field) to find the free-energy extrema for faceted clusters far from the 
critical point (e = 0.5), where even the large crystals are inherently faceted. Under 
such conditions, the free-energy surface has many local minima [28, 142] that map 
out the shape of the free-energy barrier for nucleation as a function of cluster size 
[see figure 23(a)]. It has also been reported that the effective interfacial free energy 
deduced from the barrier height using the hexagonal classical cluster model (valid 
owing to the interface sharp on the atomic scale) converges towards the free energy 
of the planar interface as the supersaturation decreases [28]. 

Toth et al. [28] have performed a similar ELE analysis in 3D to study for het- 
erophase fluctuations forming in the IM-PFC model. They have evaluated the 
properties of crystal nuclei for bcc and fee structures. It has been found that un- 
der the investigated conditions both the nucleation barrier and the driving force 
are fairly close for these structures, indicating comparable interfacial free energies 
and also Turnbull's coefficient^ for the bcc and fee structures under conditions 
(e = 0.3748), where the thermodynamic properties of the crystalline phases are 
rather close [42] . This seems to contradict recent results from MD simulations per- 
formed using the EAP method [198]. It is worth noting, however, that Turnbull's 
coefficient varies with the type of interaction potential. For example, it is ^ 0.55 
for EAP metals of fee structure [198], whereas ^ 0.36 has been deduced for the fee 
Lennard- Jones (LJ) system [199]. 

4.1.2. Heterogeneous nucleation 

The height of the nucleation barrier is often reduced by heterogeneities (walls, 
floating particles, templates, etc.), a phenomenon termed heterogeneous nucleation 



-•^ Turnbull's coefficient Ct is a reduced liquid-solid interfacial free energy defined via the relationship 
7ls = C't — w3^2/3 7 where 7is, Aiff, A^a, and Vm are the total liquid-solid interfacial free energy, the molar 

heat of fusion, the Avogadro number, and the molar volume, respectively. Ct is expected to depend only 
on the crystal structure. Recent results indicate that besides structure the interaction potential also has 
influence on its magnitude. 
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Figure 24. Heterogeneous crystal nucleation at a wall in two spatial dimensions from solving the ELE 
[140]: (a), (b) typical (non-faceted) nuclei obtained for e = 0.25, ipo = —0.341, while the ratio of the lattice 
constant of the substrate and the particle diameter is ag/cr = 1.49 and 2.0, respectively, whereas the crystal 
orientations (112) and (Oil) are parallel with the wall. The intersection of the circular and linear fits (white 
lines) to the contour line (green line) defines the contact angle, (c) Contact angle vs. ds/cr for e = 0.25 and 
ipo = —0.341. The red symbols indicate the states corresponding to panels (a) and (b). (Reproduced from 
G. I. Toth, G. Tegze, T. Pusztai, and L. Granasy, Heterogeneous crystal nucleation: the effect of lattice 
mismatch, Phys. Rev. Lett. 108 (2012), 025502, no. 2, DOI: 10. 1103/PhysRevLett. 108.025502 © 2012 by 
the American Physical Society.) 



[200]. The efficiency of the heterogeneities in instigating freezing is influenced by 
a range of microscopic properties, such as crystal structure, lattice mismatch, sur- 
face roughness, adsorption, etc., which are often condensed into the contact angle 
used in the classical theory and coarse-grained continuum models. The atomic scale 
characteristics of the substrate surface, especially the lattice mismatch, are crucial 
from the viewpoint of the highly successful free-growth-limited model of particle- 
induced freezing by Greer and co-workers [200, 201] - a model in which cylindrical 
particles, whose circular faces (of radius R) are ideally wet by the crystal, remain 
dormant during cooling until the radius of the homogeneous nuclei becomes smaller 
than R and free growth sets in. The PFC models are especially suitable to inves- 
tigate such problems as they work on the diffusive time scale [12] and can handle 
systems containing as many as 2.4 x 10^ atoms [28]. 

Along this line, Toth et al. [28] used a periodic external potential to incorporate 
a crystalline substrate into the ELE method for determining the properties of 
faceted heterogeneous nuclei. They have observed the adsorption of a monolayer 
of particles on the surface of substrate that reduced the formation energy of nuclei 
substantially and lead to a contact angle of 60° determined by the crystal structure 
[see figure 23(b)]. In a more recent ELE study, Toth etal. [140] have shown for the 
IM-PFC model in two spatial dimensions that the contact angle, the thickness of 
the crystal layer adsorbed on the substrate, and the height of the nucleation barrier 
vary non-monotonically with the lattice constant of a square-lattice substrate (see 
figure 24). They have also proven in two and three spatial dimensions that the free- 
growth-limited model of particle-induced freezing by Greer et al. [200, 201] is valid 
for larger nanoparticles and a small anisotropy of the interfacial free energy (see 
figure 25). Faceting due to either the small size of the foreign particle or a high 
anisotropy of the free energy of the liquid-solid interfacial free energy decouples 
free growth from the critical size of homogeneous nuclei. 



4.2. Pattern formation 

Owing to the overdamped diffusive dynamics most of the PFC models assume, 
diffusional instabilities that lead to fingering of the propagating crystal front are 
inherently incorporated. In the case of a single-component PFC model, diffusive dy- 
namics means that as the growing crystal (of larger particle density than the liquid) 
consumes the particles in the adjacent liquid; the only way they can be replenished 
is via long-range diffusion from the bulk liquid. Accordingly, a depletion zone forms 
ahead of the growing crystal [28, 38, 42, 202]. This resembles the behaviour of col- 
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Figure 25. Free-growth hmited bcc crystahization in 3D on a cube of sc structure as predicted 
by the ELE within the IM-PFC model [140]. Here, e = 0.25 and from left to right = 
-0.3538, -0.3516, -0.3504, -0.3489, -0.3482, and -0.3480, respectively. The linear size of the substrate 
is Ls = 16abcc5 while dbcc is the lattice constant of the stable bcc structure. Spheres centred on density 
peaks are shown, whose size increases with the height of the peak. Colour varies with peak height, inter- 
polating between red (minimum height) and light grey (maximum height). The ELE has been solved on 
a 256 X 256 x 256 grid. (Reproduced from G. I. Toth, G. Tegze, T. Pusztai, and L. Granasy, Heteroge- 
neous crystal nucleation: the effect of lattice mismatch^ Phys. Rev. Lett. 108 (2012), 025502, no. 2, DOI: 
10.1103/PhysRevLett. 108. 025502 © 2012 by the American Physical Society.) 

loidal suspensions, in which the micron-sized colloid particles move by Brownian 
motion in the carrier fluid. Relying on this similarity, the single-component PFC 
models can be considered as reasonable tools to address colloidal crystal aggrega- 
tion. One may though get rid of this type of mass-diffusion-controlled dynamics 
when driving the system strongly enough for a diffusion-controlled to diffusionless 
transition [38] . These phenomena are not necessarily present in the phenomenolog- 
ical coarse-grained PFC models, where the change in density upon crystallization 
is not always taken into account - models that might be considered, therefore, as 
a reasonable description of metallic alloys. Diffusive dynamics is furthermore ap- 
propriate on the surface of substrates, where the adsorbed particles indeed move 
by diffusion in a periodic potential. 2D PFC models relying on the interplay of the 
inter-particle forces and a periodic external potential representing the symmetries 
of the substrate can be used to capture pattern formation on such surfaces. 

4.2.1. PFC modelling of surface patterns 

A combination of computer simulations with the solution of amplitude equations 
for PFC models equipped with periodic potentials has been used to investigate 
a range of surface phenomena including incommensurate to commensurate tran- 
sitions for triangular surface layers on a square-lattice substrate [183, 184, 203], 
the deposition of monolayers on quasi-crystalline substrates [204], sliding friction 
[182, 184], and the formation of quantum dots/islands on nanomembranes [205, 206] 
(see figure 26). 

Muralidharan and Haataja [48] have extended the PFC model for describing 
stress-induced alloying of bulk-immiscible binary systems on a substrate by adding 
a potential energy term describing the substrate and a regular solution term. Fix- 
ing the model parameters to data for CoAg/Ru(0001), they demonstrated that 
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Figure 26. Growth of an island on a nanomembrane of approximately 30 atomic layers thickness at times 
(a) t = 60000, (b) t = 80000, (c) t = 400000 as predicted by the amplitude equations of reference [206]. 
Similar to figure 9, the left, middle, and right columns correspond to the results of the total particle 
density, the difference of the particle densities for the two components, and the local free-energy density, 
respectively. Also for clarity, the particle density images have been expanded by a factor of two. (Repro- 
duced with permission from K. R. Elder and Z.-F. Huang, A phase field crystal study of epitaxial island 
formation on nanomembranes, J. Phys.: Condens. Matter 22 (2010), 364103, no. 36, DOI: 10.1088/0953- 
8984/22/36/364103 © 2010 by Institute of Physics Publishing.) 




Figure 27. Surface patterns the amplitude equations predict for a decreasing dimensionless coupling uq 
between the layer and the substrate [49]. From left to right: uq = 12.1 x 10~^, 3.2 x 10~^, and 0.87 x 10~^, 
respectively. Colouring: fee domains are blue, hep domains are red, and the domain walls are green. 
(Reproduced with permission from K. R. Elder, G. Rossi, P. Kanerva, F. Sanches, S. C. Ying, E. Granato, 
C. V. Achim, and T. Ala-Nissila, Patterning of heteroepitaxial overlayers from nano to micron scales, Phys. 
Rev. Lett. 108 (2012), 226102, no. 22, DOI: 10.1103/PhysRevLett. 108. 226102 © 2012 by the American 
Physical Society.) 

the model captures experimentally observed morphologies. A similar approach has 
been proposed by Elder et al. [49] using amplitude equations that allow large-scale 
simulations for stress-induced alloying in heteroepitaxial overlayers. Quantitative 
predictions, that are in an excellent agreement with experiments, have been ob- 
tained for the stripe, honeycomb, and triangular superstructures emerging in the 
metal/metal systems, Cu on Ru(OOOl) and Cu on Pd(lll) (see figure 27). 
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Figure 28. Growth of five dendrites in the binary PFC model (the distribution of the field -0 is shown). The 
snapshots taken at 1000, 5000, 10000 and 20000 time steps are shown. The simulations have been performed 
on a 16384 x 16384 grid, using a semi-implicit spectral method [164]. (Reproduced from T. Pusztai, G. 
Tegze, G. I. Toth, L. Kornyei, G. Bansel, Z. Fan, and L. Granasy, Phase-field approach to poly crystalline 
solidification including heterogeneous and homogeneous nucleation, J. Phys.: Condens. Matter 20 (2008), 
404205, no. 40, DDI: 10.1088/0953-8984/20/40/404205 © 2008 by Institute of Physics Publishing.) 




Figure 29. Morphological transition from dendritic needle crystals to compact hexagonal shape with 
increasing driving force for crystallization [29]. Conditions/properties are as described for dendrites in 
reference [164] except that the initial total number densities are ipo = 0.009, 0.0092, 0.0094, 0.0096, and 
0.0098 (from left to right). The reduced number density ip is shown. Notice the reducing contrast of the 
images from left to right indicating an increasing solute trapping. A 8192 x 8192 grid has been used. (Re- 
produced from G. Tegze, Application of atomistic phase-field methods to complex solidification problems, 
PhD Thesis, Eotvos University, Budapest, Hungary (2009).) 



4.2.2. Pattern formation during binary solidification 

4.2.2.1. Dendritic freezing. The possibility of growing solutal dendrites has been 
first addressed within numerical IM-PFC simulations by Elder etal. [26]. Studies 
of the transformation kinetics for many particles including several dendrites have 
been performed using the same approach by Pusztai et al. [27] for system sizes 
containing about 1.6 x 10^ atoms (see figure 28). Tegze [29] has investigated the 
behaviour of solutal dendrites using binary IM-PFC simulations of similar size. 
With increasing driving force obtained by increasing the total number density, 
transitions from dendritic needle crystals to compact hexagon shape crystals have 
been observed as in the conventional PF models (see figure 29). It has also been 
shown that (i) a steady state tip velocity is attained after a time, and (ii) tip 
oscillations do not occur, i.e., from the viewpoint of side branch formation the 
dendrite tip works like a selective amplifier of the fluctuations at the tip (see figure 
30). 



4.2.2.2. Eutectic solidification. Eutectic solidification in binary IM-PFC simu- 
lations has been first observed in two spatial dimensions in the seminal paper by 
Elder etal. [26]. The formation of lamellar eutectic grains has been explored by 
Elder, Huang, and Provatas [31] using and approach based on amplitude equa- 
tions (see figure 31). For the relatively large lattice mismatch they assumed for 
the two crystalline phases (8.4% in equilibrium), spontaneous nucleation of dislo- 
cations at the lamellar interfaces has been observed - a phenomenon expected to 
modify the spacing selection mechanism predicted by earlier eutectic solidification 
theories. Larger-scale binary IM-PFC simulations relying on a numerical solution 
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(a) (b) (c) 

Figure 30. Analysis of a solutal dendrite grown in the binary IM-PFC model [29]. The dendrite arms 
are numbered clockwise from the top arm. Apparently, there are no tip radius or velocity oscillations and 
steady state growth is reached after ^ 4000 time steps. A 8192 x 8192 grid has been used and ipo = 0.0092, 
whereas other conditions as for figure 25 are present. (Reproduced from G. Tegze, Application of atomistic 
phase-field methods to complex solidification problems, PhD Thesis, Eotvos University, Budapest, Hungary 
(2009).) 




Figure 31. Time evolution of equiaxed eutectic solidification within the amplitude equation formalism 
proposed by Elder, Huang, and Provatas [31]. Panels (a), (b), and (c) correspond to dimensionless times 
30000, 60000, and 105000, respectively. From left to right the columns display the reduced total number 
density in the boxed region, the coarse-grained number density, the reduced difference of the number 
densities for the two species, and the local free-energy density. Dislocations appear as small black dots 
in the local free-energy density. (Reproduced with permission from K. R. Elder, Z.-F. Huang, and N. 
Provatas, Amplitude expansion of the binary phase- field- crystal model, Phys. Rev. E 81 (2010), 011602, 
no. 1, DOI: 10. 1103/PhysRevE.81. 011602 © 2010 by the American Physical Society.) 

of the EOMs (73) and (74) in two spatial dimensions imply owing to the diffu- 
sive dynamics of the total particle density in the binary IIM-PFC model, eutectic 
colonies may form even in binary systems [28] (see figure 32). Here, the morpho- 
logical change occurs as a result of the diffusional instability emerging from the 
diffusive EOM, many of the PFC models assume. Using the same approach in 3D, 
eutectic crystallization to the bcc structure has been reported by Toth et al. [28] 
(see figure 33). These atomistic simulations indicate a remarkable time evolution 
of the eutectic pattern after solidification. 
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(a) (b) (c) 

Figure 32. Snapshots of eutectic sohdification on the atomistic scale in the binary IM-PFC model in two 
spatial dimensions [28]: composition maps corresponding to 2x 10^, 6x 10^, and 10^ time steps are shown. 
White and black denote the two crystalline phases, while yellow stands for the liquid phase. The simulation 
has been performed on a 2048 x 1024 rectangular grid. Crystallization has been started by placing a row of 
supercritical crystalline clusters of alternating composition into the simulation window. Interestingly, the 
eutectic pattern evolves inside the solid region on a timescale comparable to the timescale of solidification. 
(Reproduced from G. I. Toth, G. Tegze, T. Pusztai, G. Toth, and L. Granasy, Polymorphism, crystal 
nucleation and growth in the phase-field crystal model in 2D and 3D, J. Phys.: Condens. Matter 22 (2010), 
364101, no. 36, DOI: 10.1088/0953-8984/22/36/364101 © 2010 by Institute of Physics Publishing.) 



•• IP 

(a) (b) (c) (d) 

Figure 33. Snapshots of eutectic solidification as predicted by the binary IM-PFC model in 3D [28]: time 
elapses from left to right. The simulation has been performed on a 450 x 300 x 300 rectangular grid. Solid- 
ification has been started by placing two touching supercritical bcc clusters of different compositions into 
the simulation window. Remarkably, the nanoscale solid-phase eutectic pattern roughens on a timescale 
comparable to the time of solidification. Brown and grey colours denote the terminal solutions of the two 
crystalline phases. Spheres of size refiecting the height of the total number density peak and coloured 
according to the local composition ip are centred to the particle density maxima. Only half of the simu- 
lation window is shown. (Reproduced from G. I. Toth, G. Tegze, T. Pusztai, G. Toth, and L. Granasy, 
Polymorphism, crystal nucleation and growth in the phase-field crystal model in 2D and 3D, J. Phys.: 
Condens. Matter 22 (2010), 364101, no. 36, DOI: 10.1088/0953-8984/22/36/364101 © 2010 by Institute 
of Physics Publishing.) 

4.3. Phenomena in the solid state 

One of the most successful areas, where the PFC models make a real difference is the 
modelling of solid state transitions including grain-boundary dynamics, melting, 
crack formation, stress-induced morphology evolution, and the modelling of the 
Kirkendall effect, to mention a few. 

4.3.1. Dislocation dynamics and grain-boundary melting 

Already the first paper on the PFC method [12, 50] has addressed grain bound- 
aries and shown that the model automatically recovers the Read-Shockley rela- 
tionship between grain-boundary energy and misorientation (see figure 34). It has 
also been shown that PFC models are ideal for modelling grain-boundary dynam- 
ics [12, 50] (see figure 34) and offers the possibility to link mechanical properties 
with the grain structure [50]. Two mechanisms of dislocation glide have been ob- 
served: for high strain rates, continuous glide is observed, while at lower strain 
rate the dislocation set into a stick-slip motion [51]. Grain-boundary melting has 
been addressed in several works [36, 37, 167]. It has been reported that disloca- 
tions in low-angle to intermediate-angle grain boundaries melt similarly until an 
angle-dependent first-order wetting transition occurs, when neighbouring melted 
regions coalesce. In the large- angle limit, the grain-boundary energy becomes in- 
creasingly uniform along its length and can no longer be interpreted in terms of 
individual dislocations (see figure 34). The difference between high- and low-angle 
boundaries appears to be reflected in the dependence of the disjoining potential 
on the width of the pre-melted layer it;: it is purely repulsive for all widths for 
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Figure 34. PFC modelling of defect and pattern formation in solids. From left to right: (a) grain-boundary 
energy vs. Read-Shockley relationship and (b) grain-boundary dynamics. (Reproduced with permission 
from K. R. Elder and M. Grant, Modeling elastic and plastic deformations in nonequilibrium processing 
using phase field crystals, Phys. Rev. E 70 (2004), 051605, no. 5, DOI: 10.1103/PhysRevB. 75. 064107 © 
2004 by the American Physical Society.) (c) Grain-boundary melting at a large-angle grain boundary. 
(Reproduced with permission from J. Mellenthin, A. Karma, and M. Plapp, Phase-field crystal study of 
grain-boundary premelting, Phys. Rev. B 78 (2008), 184110, no. 18, DOI: 10. 1103/PhysRevB. 78. 184110 © 
2008 by the American Physical Society.) 




Figure 35. Crack formation (left) and strain-induced epitaxial islands (right) in the single-component IM- 
PFG model, (a), (b) Snapshots of the energy density map taken at dimensionless times 25000 and 65000. 
(Reproduced with permission from K. R. Elder and M. Grant, Modeling elastic and plastic deformations 
in nonequilibrium processing using phase field crystals, Phys. Rev. E 70 (2004), 051605, no. 5, DOI: 
10. 1103/PhysRevB. 75. 064107 © 2004 by the American Physical Society.) (c) Gray scale image of epitaxial 
islands in an IM-PFC simulation for a 4.8% tensile film. (Reproduced with permission from Z.-F. Huang 
and K. R. Elder, Mesoscopic and microscopic modeling of island formation in strained film epitaxy, Phys. 
Rev. Lett. 101 (2008), 158701, no. 15, DOI: 10. 1103/PhysRevLett. 101. 158701 © 2008 by the American 
Physical Society.) 

misorientations larger than a critical angle, however, it switches from repulsive at 
small w to attractive for large w [37]. 

4.3.2. Crack formation and propagation 

Elder and Grant [50] have demonstrated that the IM-PFC model can be used 
to model crack propagation. A small notch cut out of a defect-free crystal placed 
under 10% strain in the vertical direction and filled with coexisting liquid has been 
used as a nucleation cite for crack propagation. Snapshots of crack development 
are shown in figure 35. 

4.3.3. Strain-induced morphologies 

Huang and Elder [52] have studied strain-induced film instability and island 
formation using numerical IM-PFC simulations and amplitude equations (see fig- 
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Figure 36. Density trapping as predicted by the single-component IM-PFC model [38, 39]. (a) Coarse- 
grained particle densities ip for the liquid and solid phases at the growth front as a function of growth 
velocity v. (b) Effective partition coefficient k defined using the liquidus and solidus densities vs. growth 
velocity. For comparison, fits of the models by Aziz [207] and Jackson etal. [208] are also displayed, 
(c) Comparison of the interface thickness d and the diffusion length di) as a function of growth velocity. 
(Reproduced from G. Tegze, L. Granasy, G. I. Toth, J. F. Douglas, and T. Pusztai, Tuning the structure of 
nonequilibrium soft materials by varying the thermodynamic driving force for crystal ordering, Soft Matter 
7 (2011), 1789-1799, no. 5, DOI: 10.1039/c0sm00944j © 2011 by Royal Society of Chemistry Publishing.) 

ure 34). They have identified a linear regime for the island wave number scaling 
and recovered the continuum ATG instability in the weak strain limit. The ATG 
instability has been studied in two spatial dimensions by Spatschek and Karma 
[54] using a different amplitude equations approach. Qualitatively similar surface 
roughening has been reported by Tegze et al. [42] for heteroepitaxial body-centred 
tetragonal (bet) films grown on sc crystalline substrates of tuned lattice constant 
- a phenomenon interpreted in terms of the Mullins-Sekerka/ATG instability. 

4.3.4- Kirkendall effect 

Elder, Thornton, and Hoyt [55] have used a simple extension of the binary IM- 
PFC model incorporating unequal atomic mobilities to investigate different as- 
pects of the Kirkendall effect. They have shown that the model indeed captures 
such phenomena as crystal (centre-of-mass) motion, pore formation via vacancy 
supersaturation, and enhanced vacancy concentration near grain boundaries. 

4.3.5. Density /solute trapping 

In recent works by Tegze etal. [38, 39], it has been reported for the IM-PFC 
model (of diffusive dynamics) that at large e (= 0.5) and high driving force a 
transition from diffusion-controlled to diffusionless solidification can be observed, 
during which the interface thickness increases, whereas the density difference be- 
tween the crystal and the liquid decreases drastically (see figure 36). This "density 
trapping" phenomenon is analogous to solute trapping observed in rapid solidifica- 
tion of alloys (where due to a lack of time for partitioning, solids of nonequilibrium 
compositions form) and can be fitted reasonably well using the models of Aziz [207] 
and Jackson etal. [208]. In a very recent work, Humadi, Hoyt, and Provatas [45] 
have investigated solute trapping in the binary MPFC model. In agreement with 
the findings for density trapping, they have found that pure diffusive dynamics 
leads to a velocity-dependent partition coefficient that approaches unity for large 
velocities - consistently with the model of Aziz and Kaplan [209]. In contrast, the 
wavelike dynamics, the second-order time derivatives of the MPFC-type EOMs 
realize, leads to a solute trapping behaviour similar to the predictions of Galenko 
etal. [210]. 

4.3.6. Vacancy / atom transport in the VPFC model 

The VPFC model by Chan, Goldenfeld, and Dantzig [56] is one of the most ex- 
citing extensions of the original PFC approach. The extra term added to the free 
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Figure 37. VPFC modelling of fluid and crystalline states of different particle densities. The number of 
atoms increases from left to right and from top to bottom. (Reproduced with permission from P. Y. Chan, 
N. Goldenfeld, and J. Dantzig, Molecular dynamics on diffusive time scales from the phase- field- crystal 
equation, Phys. Rev. E 79 (2009), 035701R, no. 3, DOI: 10.1103/PhysRevE.79.035701 © 2009 by the 
American Physical Society.) 

energy makes particle density non-negative and allows for the formation of individ- 
ual density peaks ("atoms" forming the fluid) and vacancies in the crystal. This, 
combined with the MPFC EOM (64), that considers inertia and damping, makes it 
a kind of MD-like approach working on a still far longer time scale than the usual 
MD simulations. Accordingly, one can obtain configurations that look like snap- 
shots of the fluid state (see figure 37) and may evaluate the structure factor for the 
fluid state, which is evidently impossible for the original PFC model. (Apparently, 
similar images can be obtained in the IM-PFC model as a transient state during 
solidification [38], however, with different dynamics owing to the differences in free 
energy and EOM.) 

A comparison with another recent development, termed the diffusive molecular 
dynamics (DMD) technique, by Li etal. [211] would be very interesting. The latter 
approach works on the diffusive time scale too, while maintaining atomic resolution, 
by coarse-graining over atomic vibrations and evolving a smooth site-probability 
representation. 

5. Phase- field-crystal modelling in soft matter physics 
5.1. Applications to colloids 

In this section, we review results obtained using different PFC models relying on 
overdamped conservative dynamics - a reasonable approximation for colloidal crys- 
tal aggregation. We concentrate on three major areas: crystal nucleation, pattern 
formation in free growth, and pattern formation in the presence of external poten- 
tials. 

As mentioned previously, using of the EOM for simulating crystallization is 
not without difficulties. In the DDFT-type models, the system cannot leave a 
metastable state (e.g., the homogeneous initial fluid) unless Langevin noise repre- 
senting thermal fluctuations is added to the EOM. This raises, however, essential 
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questions: considering the number density an ensemble-averaged quantity, all the 
fluctuations are (in principle) incorporated into the free-energy functional. Adding 
noise to the EOM, a part of the fluctuations might be counted twice [95, 192]. 
If in turn the number density is viewed as being coarse-grained in time, there 
is phenomenological motivation to add a noise term to the EOM [193]. The latter 
approach is appealing in several ways: crystal nucleation is feasible from a homoge- 
neous state and capillary waves appear at the crystal-liquid interface. To investigate 
how nucleation and growth happen on the atomistic level, a conserved noise term is 
usually incorporated into the EOM [see equations (60)-(63)]. To overcome some dif- 
ficulties occurring when discretising the noise [194, 195], coloured noise obtained 
by filtering out the unphysical short wavelengths smaller than the inter-particle 
distance is often used (this removes both the ultraviolet catastrophe expected in 
3D [212] and the associated dependence of the results on spatial resolution). The 
majority of the studies we review below follows this approach. 

5.1.1. Nucleation in colloidal crystal aggregation 
5.1.1.1. Homogeneous nucleation. 

The effect of noise. A systematic study of the effect of the noise strength on 
the grain size distribution performed in two spatial dimensions by Hubert et al. [34] 
for the original IM-PFC model implies that grain size decreases with increasing 
noise amplitude, resulting in both a smaller average grain size and a reduced max- 
imum grain size. They have distinguished two regimes regarding the cluster size 
distribution: for small noise amplitudes a bimodal cluster size distribution is ob- 
served, whereas for large noise amplitudes a monotonically decreasing distribution 
is reported. 

Phase selection in 2D and 3D. Mounting evidence indicates that the classical 
picture of crystal nucleation, which considers heterophase fluctuations of only the 
stable phase, is oversimplified. Early analysis by Alexander and McTague suggests 
a preference for bcc freezing in simple liquids [213]. Atomistic simulations for the 
LJ system have verified that small heterophase fluctuations have the metastable 
bcc structure, and even larger clusters of the stable fee structure have a bcc in- 
terface layer [214], while the ratio of the two phases can be tuned by changing 
the pressure [215]. Composite bcc-fcc nuclei have also been predicted by contin- 
uum models [216]. Two-stage nucleation has been reported in systems that have 
a metastable critical point in the undercooled liquid (including solutions of glob- 
ular proteins [217]); the appearance of the crystalline phase is assisted by dense 
liquid droplets, whose formation precedes and helps crystal nucleation [218]. Re- 
cent studies indicate a similar behaviour in simple liquids such as the LJ [219] or 
hard-sphere (HS) [220] fluids, where a dense liquid or amorphous precursor assists 
crystal nucleation. Analogous behaviour has been reported for colloidal systems 
in two spatial dimensions [221] and 3D [222]. These findings imply that the nu- 
cleation precursors are fairly common. Since the IM-PFC model has bcc, fee, and 
hep stability domains [28], the appearance of an amorphous phase and two-step 
nucleation has also been reported [46], and the 2M-PFC model incorporates the 
IM-PFC model [132]. This class of the dynamic PFC models is especially suitable 
for investigating phase selection during freezing of undercooled liquids. 

In two spatial dimensions, it has been shown within the framework of the IM- 
PFC model that at relatively small supersaturations direct crystal nucleation takes 
place. Increasing the thermodynamic driving force, first copious crystal nucleation 
is observed, and at higher driving forces an amorphous precursor precedes crys- 
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(f) (g) (h) (i) (j) 

Figure 38. Snapshots of early and late stages of isothermal solidification in IM-PFC quenching 
simulations performed in two spatial dimensions with initial reduced particle densities of ipo = 
—0.55, —0.50, —0.45, —0.40 and —0.35 [33]. (a)-(e): Early stage: the respective reduced times are r/Ar = 
10000,3000,1500,1000, and 700. (f)-(j): Late stage: the same areas are shown at reduced time r/Ar = 
60000. Reduced particle density maps in 418 x 418 sized fractions of 2048 x 2048 sized simulations are 
shown. Other simulation parameters were e = 0.75 and a = 0.1 (noise strength). (Reproduced from L. 
Granasy, G. Tegze, G. I. Toth, and T. Pusztai, Phase-field crystal modelling of crystal nucleation, het- 
eroepitaxy and patterning, Philos. Mag. 91 (2011), 123-149, no. 1, DOT: 10.1080/14786435.2010.487476 © 
2011 by Taylor & Francis.) 




Figure 39. Structural properties evolving after quenching in IM-PFC simulations [33]: (a) pair-correlation 
function g(r) for the early-stage solidification structures shown in figures 20(b)-(e). (b) Time evolution of 
the bond-order correlation function geir) for ipo = —0.4 on log-log scale, geir) is shown at r/Ar = 
1000,4000, 16000, and 64000. For comparison, the upper envelope expected for the hexatic phase and the 
result for a single crystal are also shown. These curves describe an amorphous to polycrystalline transition 
[see figures 20(d) and 20 (i)]. Notice that the upper envelope of the gQ(r) curves decay faster than expected 
for the hexatic phase. (Reproduced from L. Granasy, G. Tegze, G. I. Toth, and T. Pusztai, Phase-field 
crystal modelling of crystal nucleation, heteroepitaxy and patterning, Philos. Mag. 91 (2011), 123-149, no. 
1, DOT: 10.1080/14786435.2010.487476 © 2011 by Taylor & Francis.) 

talline nucleation [33] [see figures 38 and 39(a)]. Similarly to quenching experiments 
for two spatial dimensions colloidal systems [223], no hexatic phase is observed in 
the IM-PFC quenching simulations [33, 38] [as demonstrated by the form of the 
radial decay of the bond-order correlation function [33], see figure 39(b)]. 

In 3D, a systematic dynamic study of the 1M/2M-PFC models by Toth et al. [35] 
shows that in these systems the first appearing solid is amorphous, which promotes 
the nucleation of bcc crystals (see figure 40) but suppresses the appearance of the 
fee and hep phases. The amorphous phase appears to coexist with the liquid in- 
dicating a first-order phase transition between these phases in agreement with the 
observed nucleation of the amorphous state. Independent ELE studies determining 
the height of the nucleation barrier have confirmed that density and structural 
changes take place on different times scales [35]. This finding suggests that the two 
time scales are probably present independently of the type of dynamics assumed. 
These findings have been associated with features of the effective interaction po- 
tential deduced from the amorphous structure using Schommer's iterative method 
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Figure 40. Two-step nucleation in the IM-PFC model at -00 = -0.1667 and e = 0.25 [35]. Four pairs 
of panels are shown, where t indicates the time elapsed, N is the total number of particles, and qi is 
the bond orientational order parameter with index i. Left: snapshots of the density distribution taken 
at the dimensionless times r = 57.74f. Spheres of the diameter of the inter-particle distance centred on 
density peaks higher than a threshold (= 0.15) are shown. They are coloured red if q4 £ [0.02; 0.07] and 
qq G [0.48; 0.52] (bcc-like) and light grey otherwise. Right: population distribution of qq (histogram painted 
similarly) and the time dependence of the fraction X of bcc-like neighbourhoods (solid line). (Reproduced 
from G. I. Toth, T. Pusztai, G. Tegze, G. Toth, and L. Granasy, Amorphous nucleation precursor in highly 
nonequilibrium fluids, Phys. Rev. Lett. 107 (2011), 175702, no. 17, DOL 10. 1103/PhysRevLett. 107. 175702 
© 2011 by the American Physical Society.) 

[224] that shows a maximum at ro\/2, where tq is the radius corresponding to the 
main minimum of the potential. Such a maximum in the interaction potential is ex- 
pected to suppress crystalhzation to the close-packed structures fee and hep [225], 
whereas the multiple minima also found are expected to lead to coexisting disor- 
dered structures [226]. By combining the results available for various potentials 
(LJ [219], HS [220], and the PFC potentials [28, 35]), it appears that a repulsive 
core suffices for the appearance of a disordered precursor, whereas the peak at 
ro\/2 correlates with the observed suppression of fee and hep structures, while the 
coexistence of the liquid and amorphous phases seen here can be associated with 
multiple minima of the interaction potential. 

3D studies, performed for bee crystal nucleation in molten pure Fe in the frame- 
work of the EOF-PFC model [28, 33], lead to similar results, however, still with 
diffusive dynamics. In these simulations, the initial density of the liquid has been 
increased until the solidification started - a procedure that has lead to an extreme 
compression owing to the small size and short time accessible for the simulations. 
While this raises some doubts regarding the validity of the applied approximations, 
the behaviour observed for the EOF-PFC Fe is fully consistent with the results ob- 
tained for the IM-PFC model: with increasing driving force first an amorphous 
precursor nucleates and the bee phase appears inside these amorphous regions 
[28, 33]. At higher driving forces the amorphous precursor appears nearly homo- 
geneously in space and the bee phase nucleates into it later. Apparently, direct 
nucleation of the bee phase from the liquid phase requires a longer time than via 
the amorphous precursor, suggesting that the appearance of the bee phase is as- 
sisted by the presence of the amorphous phase and in line with recent predictions 
by DFT [219] and atomistic simulations [220]. Remarkably, the interaction poten- 
tial evaluated for Fe from the pair-correlation function of the amorphous structure 
is oscillatory and is qualitatively similar to the ones evaluated from experimental 
liquid structures [227]. 

5.1.1.2. Heterogeneous nucleation. Prieler et al. [137] have explored crystal nu- 
cleation on an unstructured hard wall in an anisotropic version of the IM-PFC 
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Figure 41. Heterogeneous nuclei formed on a hard wall in the APFC model proposed in reference [137] and 
the dependence of the left and right side contact angle (71 and 72, respectively) on the crystal orientation. 
(Reproduced from R. Prieler, J. Hubert, D. Li, B. Verleye, R. Haberkern, and H. Emmerich, An anisotropic 
phase-field crystal model for heterogeneous nucleation of ellipsoidal colloids, J. Phys.: Condens. Matter 21 
(2009), 464110, no. 46, DOI: 10.1088/0953-8984/21/46/464110 © 2009 by Institute of Physics Publishing.) 



model, in which the particles are assumed to have an ellipsoidal shape. In par- 
ticular, they have investigated how the contact angle depends on the orientation 
of the ellipsoids and the strength of the wall potential (see figure 41). A complex 
behaviour has been observed for the orientational dependence, while increasing the 
strength of the wall potential reduced the contact angle. 

Granasy et al. [136] have studied crystal nucleation in an rectangular corner of 
structured and unstructured substrates within the IM-PFC model in two spatial 
dimensions. Despite expectations based on the classical theory of heterogeneous nu- 
cleation and conventional PF simulations [228], which predict that a corner should 
be a preferred nucleation site, in the atomistic approach such a corner is not a 
preferable site for the nucleation of the triangular crystal structure (see figure 42) 
owing to the misfit of the triangular crystal structure with a rectangular corner. 
Crystals of different orientation nucleate on the two substrate surfaces, which in- 
evitably leads to the formation of a grain boundary starting from the corner when 
the two orientations meet. The energy cost of forming the grain boundary makes 
the rectangular corner an unfavoured place for nucleation. In contrast, a 60° corner 
helps the nucleation of the triangular phase. 

5.1.2. Pattern formation in colloidal crystal aggregation 

5.1.2.1. Colloid patterns in two dimensions. Choosing a large value for the 
parameter e where the liquid-solid interface is faceted, Tegze etal. [38, 39] have 
investigated solidification morpholog les as a function of the thermodynamic driv- 
ing force. It has been found that the diffusion-controlled growth mode observed at 
low driving forces and characterized by faceted interfaces changes to a diffusion- 
less growth mode of a diffuse liquid-solid interface that produces a crystal, whose 
density is comparable to the density of the liquid due to quenched-in vacancies 
(see section 4.3.5). These two modes have already been observed experimentally in 
colloidal systems [229] . It has been shown that the two modes can coexist and lead 
to a new branching mechanism that differs from the usual diffusional instability 
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(a) (b) (c) 

Figure 42. Heterogeneous nucleation in rectangular inner corners of the IM-PFC model in two spatial 
dimensions [33]. (a) Nucleation on (01) surfaces of a square lattice (ratio of lattice constant of substrate 
to inter-particle distance ao/a ^ 1.39). (b) Nucleation on (11) surfaces of a square lattice, (c) Nucleation 
on an unstructured substrate. Notice the frustration at the corner and the formation of a grain boundary 
starting from the corner at later stages. (Reproduced from L. Granasy, G. Tegze, G. I. Toth, and T. 
Pusztai, Phase-field crystal modelling of crystal nucleation, heteroepitaxy and patterning, Philos. Mag. 91 
(2011), 123-149, no. 1, DOT: 10.1080/14786435.2010.487476 © 2011 by Taylor & Francis.) 




Figure 43. Single crystal growth morphologies (a)-(d) in the IM-PFC model [38] (top) and experiment 
(bottom) (e)-(h): 2D colloid crystals by Skjeltorp. (Reproduced with permission from A. T. Skjeltorp, 
Visualization and characterization of colloidal growth from ramified to faceted structures, Phys. Rev. 
Lett. 58 (1987), 1444-1447, no. 14, DOT: 10.1103/PhysRevLett.58.1444 © 1987 by the American Physical 
Society.) The driving force increases from left to right. In the case of the simulations, the coarse-grained 
particle density map is shown. The fractal dimensions of the single crystal aggregates evaluated from the 
slope of the plot \og{N) vs. \og{Rg) (N is the number of particles in the cluster and Rg is its radius of 
gyration) are: (a) /a = 2.012 ±0.3%, (b) 1.967 ±0.3%, (c) 1.536 ±0.9%, (d) 1.895 ±0.3%. The fast growth 
mode is recognizable via the lack of a (dark) depletion zone at the interface, whose presence is indicative 
to the slow mode. A 2048 x 2048 rectangular grid corresponding to ~ 13000 particles, or 118 /xm x 118 /im 
(assuming 1.1 /im particles) has been used - a size comparable to that shown by the experimental images. 
(Reproduced from G. Tegze, L. Granasy, G. I. Toth, J. F. Douglas, and T. Pusztai, Tuning the structure of 
nonequilibrium soft materials by varying the thermodynamic driving force for crystal ordering, Soft Matter 
7 (2011), 1789-1799, no. 5, DOT: 10.1039/c0sm00944j © 2011 by Royal Society of Chemistry Publishing.) 

driven branching by which dendritic structures form. This new mechanism explains 
the fractal-hke and porous growth morphologies [230] observed in 2D colloidal sys- 
tems (see figure 43) and may be relevant for the diffusion-controlled to diffusionless 
transition of crystallization in organic glasses. 



5.1.2.2. Colloid patterns in three dimensions. Toth et al. [28] have demon- 
strated first that owing to the conservative dynamics, the EOM of the IM-PFC 
model realizes, dendritic growth forms of bcc and fee structure evolve in the single- 
component theory. Tegze [29] and Granasy etal. [231] have shown by simulations 
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Figure 44. 3D crystal growth morphologies grown from a bcc seed in the single-component IM-PFC model 
at e = 0.3748 in a system containing about 3 x 10^ colloidal particles [29]. The initial fluid density decreases 
as (a) ipo = -0.015, (b) -0.0175, (c) -0.01875, (d) -0.02, (e) -0.02062, (f) -0.0225, (g) -0.025, (h) -0.03, 
(i) —0.0325. The simulations have been performed on a 1024 x 1024 x 1024 grid. Assuming 1 fim diameter 
for the particles, the linear size of the simulation box is ~ 0.16 mm - comparable to the smaller colloidal 
dendrites seen in microgravity experiments [233]. (Reproduced from G. Tegze, Application of atomistic 
phase-field methods to complex solidification problems, PhD Thesis, Eotvos University, Budapest, Hungary 
(2009).) 

containing ^3x10^ particles that due to a kinetic roughening of the crystal-hquid 
interface that leads to interface broadening, a transition can be seen from faceted 
dendrites to compact rounded crystals (see figure 44) - a phenomenon reported 
earlier in experiments for dendritic growth of NH4Br crystals [232]. Notice that 
such a kinetic effect cannot be easily incorporated into conventional PF models. 
Remarkably, as pointed out in reference [28] , assuming a micrometer diameter for 
the "atoms" , these dendritic structures are comparable in size to those formed in 
colloid experiments in microgravity [233]. This is a unique situation indeed: an 
"atomistic" theory works here on the size scale of experimental dendrites. 

In a recent work. Tang et al. [30] have performed a geometric analysis of bcc 
and fee dendrites grown in the respective stability domains of the IM-PFC model, 
and evaluated dynamic exponents characterizing dendritic growth in the (100), 
(110), and (111) directions. They associate the relatively large values obtained 
for the stability constant from the geometry of the dendrite tip with the faceted 
morphology of the crystals. 
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Figure 45. (a) Single and multiple occupation of a chemically patterned periodic substrate by colloidal 
particles as a function of increasing patch size in the experiments. (Reproduced with permission from I. 
Lee, H. Zheng, M. F. Rubner, and P. T. Hammond, Controlled cluster size in patterned particle arrays 
via directed adsorption on confined surfaces, Adv. Mater. 14 (2002), 572-577, no. 8, DOI: 10.1002/1521- 
4095(20020418)14:8<572::AID-ADMA572>3.0.CO;2-B © 2002 by Wiley.) (b) IM-PFC simulations [33] 
with increasing diameter of circular attractive potential wells. Reduced particle density maps are shown. 
The ratio of the potential well diameters relative to the single occupation case has been 1, 1.25, 1.5, 2, 2.13, 
and 2.5. (Reproduced from L. Granasy, G. Tegze, G. I. Toth, and T. Pusztai, Phase- field crystal mod- 
elling of crystal nucleation, heteroepitaxy and patterning, Philos. Mag. 91 (2011), 123-149, no. 1, DOI: 
10.1080/14786435.2010.487476 © 2011 by Taylor & Francis.) 

5.1.3. Colloid patterning 

Colloid patterning under the influence of periodic substrates can be realized via 
creating patches that are chemically attractive to the colloidal particles [234]. De- 
pending on the size of the patches single, double, triple, etc., occupations of the 
patches are possible (see figure 45), whereas the distance of the patches may lead 
to the formation of various ordered patterns, as predicted by Langevin simulations, 
in which the patterned substrate is represented by appropriate periodic potentials 
[235]. Granasy et al. [33] has employed a IM-PFC model supplemented with a peri- 
odic potential of circular potential wells arranged on a square lattice, to reproduce 
the patterns seen in the experiments (see figure 45). Another problem, exempli- 
fying the abilities of PFC simulations in modelling colloid patterning, is colloidal 
self-assembly under the effect of capillary-immersion forces acting on the colloid 
particles in thin liquid layers due to capillarity and a periodically varying depth 
of the liquid layer due to a wavy substrate surface. Experiments of this kind have 
been used to produce single and double particle chains [236] and the otherwise un- 
favourable square-lattice structure [237]. The capillary-immersion forces can often 
be well represented by a potential of the form U — ui cos(A:x), where ui is a con- 
stant, k = 27r/A, and A the wavelength of the periodic potential. Setting A = cr/v^, 
where a is the interparticle distance, and varying the orientation of the grooves 
relative to the crystallization (drying) front, patterns seen in the experiments [237] 
are observed to form in the IM-PFC model: for grooves parallel to the front, a 
frustrated triangular structure of randomly alternating double and triple layers 
appears. For grooves perpendicular to the front, the particles align themselves on 
a square lattice with the (11) orientation lying in the interface, while for a 7r/4 
declination of the grooves the same structure forms, however, now with the (10) 
face lying in the front. Using larger wavelengths for the potential and adding a 
weak transversal modulation, while starting from a homogeneous initial particle 
density, nucleation and growth of wavy single and double chains resembling closely 
to the experiments [236] are seen [33] (see figure 46). 

Epitaxial growth on the (100) surface of a sc substrate has been investigated in 
3D using the IM-PFC model by Tegze etal. [42]. The lattice constant of the 
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Figure 46. Patterning in experiment vs. IM-PFC simulation: (a) single and double particle chains evolving 
in experiment due to capillary-immersion forces on the surface of a rippled substrate. (Reproduced with 
permission from A. Mathur, A.-D. Brown, and J. Erlebacher, Self-ordering of colloidal particles in shallow 
nanoscale surface corrugations, Langmuir 22 (2006), 582-589, no. 2, DOI: 10.1021/la0520379 © 2006 by 
the American Chemical Society.) (b) The particle chains forming in the IM-PFC simulation performed with 
a tilted and wavy version of the potential described in the text [33]. Only a fraction of the reduced particle 
density map is shown. (Reproduced from L. Granasy, G. Tegze, G. I. Toth, and T. Pusztai, Phase-field 
crystal modelling of crystal nucleation, heteroepitaxy and patterning, Philos. Mag. 91 (2011), 123-149, no. 
1, DOI: 10.1080/14786435.2010.487476 © 2011 by Taylor & Francis.) 




Figure 47. In-plane snapshot of crystalline aggregates grown on 5 x 5 square-lattice templates of as/crf^c = 
1.0, 1.1547, and 1.56 in 3D as predicted by IM-PFC simulations [240]. Here, afcc = 1-056. A 256 x 256 x 128 
grid has been used. The visualization is as in figure 25. 



substrate has been varied in a range that incorporates the interatomic distance 
of the bulk fee structure and the lattice constant of the bulk bcc phase, where 
the (100) face of the sc structure is commensurable with the (100) faces of the 
bulk fee and bcc structures, respectively. A bet structure has grown, whose axial 
ratio c/a varies continuously with the lattice constant of the substrate, where c 
and a are the lattice constants of the bet structure perpendicular and parallel to 
the surface of the substrate, respectively. At the matching values of as, fee and 
bcc structures have been observed respectively, as observed in colloid patterning 
experiments [238]. Analogous results have been obtained for the (100) face of an 
fee substrate using IM-PFC simulations, however, for large lattice mismatch amor- 
phous phase mediated bcc nucleation has been seen [140]. 

Optical tweezers are used widely to realize 2D periodic templates for influenc- 
ing colloidal crystal aggregation in 3D [239]. Such templates, depending on the 
mismatch to the crystalline structure evolving, may instigate the formation of 
single-crystal or polycrystalline structures [240]. Growth textures, obtained when 
supplementing the IM-PFC model with a 5 x 5 flat square-lattice template (realized 
by a periodic potential term), show remarkable resemblance to the experiments (see 
figure 47) [241]. 
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5.2. Applications to polymers 

In the field of soft matter simulations the study of polymers plays an important 
role. Understanding the relation between the statistical physics on the molecular 
scale and the resulting material properties is of great relevance with respect to 
principal questions, chemical and pharmaceutical applications, as well as the study 
of biopolymers in living cells. The investigation of polymer systems and the devel- 
opment of new polymer materials and composites with specific properties depend 
essentially on suitable numerical and simulation techniques. One of the big chal- 
lenges is closing the gap between the atomistic and the microscale. Furthermore, 
in many polymer systems structures develop rather slowly and on very different 
length scales. Since particle-based simulations are strongly restricted, especially 
regarding the time scale, there is the demand for continuum models such as PF 
and PFC theories. 

Polymer systems are studied in various systems, in the form of polymer solutions 
and melts, in hydrogels, in glasses, and in polymer crystals. A principle problem 
of continuum models for polymers is the connectivity of the chains, which leads to 
long-ranging interactions between spatially separated atoms of the same polymer 
chain and, due to entanglements, also between atoms of different chains. For an 
implementation of these effects into continuum theories, appropriate simplifications 
are necessary. 

One major continuum method for polymer systems is the self-consistent field 
(SCF) theory. Starting from an equilibrium theory, it has been extended in order 
to study dynamic systems. It is therefore interesting to compare the SCF theory 
with the PF method, which, during the last years, has increasingly been applied on 
polymer systems, and the PFC method, which is suited for studying pronounced 
periodic structures in polymers. In this context, we will also elucidate the relation 
between the PF and the PFC method, especially for the case of polymer systems. 

In many applications of polymers, the structure and dynamics of the system are 
sensitive to various parameters like temperature, flow velocity, velocity gradients, 
gravity or stresses. In polyelectrolytes, also charge densities, ion concentrations, 
the pH value, and external fields need to be considered. Several of these aspects 
have been taken into account in PF calculations, as we will demonstrate for some 
examples. We begin with a description of the SCF theory, after which we discuss 
the PF and the PFC approaches to polymer systems. 

5.2.1. Self- consistent field theory for polymer systems 

In general, the SCF theory is a molecular-based, approximate description of 
many-body systems. For polymer blends, it has originally been developed to de- 
scribe equilibrium systems [242, 243]. Later it was extended to dynamic self- 
consistent field (DSCF) theory, which allows investigating dynamic processes (for 
an review see reference [244]). 

The SCF method is used to study fluctuations, pattern formation, and segrega- 
tion in various kinds of polymer mixtures (see references [244-247] and references 
therein). In many works, the method is applied on a mixture of different types 
of homopolymers or block copolymers built of two types of monomers A and B. 
In most cases, it is assumed that the polymer chains have a Gaussian distribu- 
tion. This means that the polymer path is distributed like a random walk and 
corresponds to the situation where the persistence length is much smaller than the 
polymer length. It should be mentioned that one can also use semiflexible, worm- 
like polymers instead [245]. Assuming an incompressible polymer melt, the total 
polymer density ^po is fixed. The canonical partition function can be written as 
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[245] 



/ Jj5r^P(r^)exp 





5((^A + ^B-^o) (119) 



Here, / 5viP{Yi) is a path integral of a polymer weighted by the Gaussian dis- 
tribution P{ri) and J^hf is the Flory-Huggins free energy for polymers, which 
favours adjacent monomers of the same type: A-A and B-B. The delta function 
generates the incompressibility. The problematic parts of Z in equation (119) are 
the polymer-polymer interactions (the steric and the Flory-Huggins interactions). 
They are avoided by rewriting >Z in a form in which single polymers are subject 
to auxiliary fields that represent the effect of the other polymers. For the system 
with monomers A and B, two auxiliary fields V and W are introduced, where V is 
conjugated to the total density (fA + ^B^ while W is conjugated to the composition 
^A — ^B' Auxiliary fields are an essential concept of all SCF theories. The canonical 
partition function Z can now be written as 



where the Hamiltonian T-L(V^W) is determined by the partition functions Zgp of 
single polymers in the presence of the auxiliary fields. The partition functions can 
be calculated with the help of propagator functions. The propagator functions are 
the solution of diffusion equations [243], which depend on V and W. Equation 
(120) is still equivalent to equation (119), but the expression is much simpler. 



5.2.1.1. Equilibrium SCF theory. In principle, one can use equation (120) 
to extract equilibrium values of the average concentration {(fA) (and thus of 
(^b) = ^0 — (^a)) and fluctuation terms. In fact, in some cases, this has been done 
numerically with the help of "field MC" simulations [248] or complex Langevin 
simulations [249]. However, for most systems the numerical cost for treating ex- 
pression (119) with fully fluctuating auxiliary fields V and W is extremely high. 
The problem can be simplified with the help of a saddle-point approximation: in- 
stead of integrating over the whole V and W space, one determines the minimum 
of the Hamiltonian T-L with respect to V and VF, 



which provides saddle-point fields V* and for given concentration fields cpA and 
(fB- In turn, the concentration fields can be calculated from the partition function 
by solving the diffusion equation that depends on V and VF*. Thus, the problem 
can be solved self-consistently by iterating the calculation of concentration and 
auxiliary fields [250]. This is the basic principle of the SCF theory. In this form, 
it enables calculation of equilibrium properties of concentration fields, including 
microphase separation and global phase separation for polymer blends including 
homopolymers and block copolymers. A more precise description is given by the 
external potential theory, which allows fluctuations of the composition field W 




(120) 



dn{v*,w*) _ dnjv*, w*) 

dV ~ dW 



= 



(121) 



[251]. 
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5.2.1.2. Dynamic SCF theory. The dynamic self-consistent field theory (DSCF) 
is an extension of the SCF theory that allows to study nonequilibrium polymer 
blends. For an incompressible fluid, the saddle point values and M^* are deter- 
mined by (pA so that we can define the energy functional 

n^A\=n{V*{ipA),W*{^A)) . (122) 
The dynamics of (pA is described by a Langevin equation 

d(pA{r,t) 



dt 



= -Vr-Jv,(r,t) + C(r,i) (123) 



with a flux density J(^(r, t) and a Gaussian random noise term ({r^t). The flux 
density is given by 



with an Onsager coefficient A^(r^ r'), which is in general nonlocal, due to the forces 
that propagate along the polymer backbone. Nonlocal effects can be considered 
with the Rouse model [252] or the reptation model [72]. In the simplest case, 
the interactions due to the chain connectivity are neglected and a local Onsager 
coefficient A(^(r,r^) = k (fA{'^^t)(pB{^^t)6{r — r') with a constant k is used [26, 85], 
leading to 

^'^^ = -kVr (fAifO - fA)Vr-l^) + C • (125) 



dt \ 5 if A 

Notice that still includes the auxiliary fields and each propagation in time re- 
quires the iterative calculation of these fields. In the random phase approximation, 
which is valid for suitably small fluctuations, one can derive an approximate ex- 
pression for J^, which has the form of a Ginzburg-Landau free energy [245], 

JbL = jdv (/fh((^a) + KVrV^)') , (126) 

with the Flory-Higgins free-energy density /fh(^a) and the parameter h. 

In analogy to DSCF, a dynamic version of the external potential (DEP) method 
can be applied. The DSCF method and the DEP theory have been used to study 
spinodal decomposition. A quantitative comparison with Monte-Carlo (MC) sim- 
ulations shows a distinctly higher accuracy of the DEP method [253]. The DSCF 
has been combined with the NS equation and a convection term has been added 
to the Langevin equation [254-256]. 

5.2.2. Phase-field methods for polymer systems 

5.2.2.1. Conventional phase-field models. In recent years, there has been a 
growing number of articles that apply the PF method to polymer systems, i.e., 
to systems in which at least one phase is a polymer melt or a polymer solution 
(see, for example, references [257-261]). Here, the PF method is used to describe 
the nonequilibrium dynamics of phase boundaries. The method is based on a free- 
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energy functional of a phase variable '^(r), which has the form [24] 

T = jdr ifi^P) + /inh(Vr^)) , (127) 

where / is the free-energy density of the bulk - usually an asymmetric double- well 
potential. The second part considers spatial inhomogeneities. A simple example of 
the free energy in equation (127) is the standard Ginzburg-Landau free energy 

J-= jdv(a{^^-^lf + h{V,^f) (128) 

with coefficients a and b. If ip is conserved, then the dynamics is described by a 
Cahn-HilUard (CH) equation [262]: 



^ = -Vr-J . (129) 



Here, the flux density is given by 



J = _MVr^ (130) 



with a mobility M, which in general may depend on ijj. 



5.2.2.2. Polymer liquids. In a fluid polymer blend, ^ is often taken as the con- 
centration of one component and the bulk free-energy density /(r) is the Margules 
activity from the Flory-Huggins theory [263, 264]. In many cases, a Gaussian noise 
term ("(r, t) is added to the CH equation so that the equation has a Langevin form 

^ = -Vr-J + C = Vr- (^Vr^^ + C • (131) 

Notice the similarity to equation (123) from SCF theory. The expressions differ in 
the free-energy functional J^, which includes the auxiliary fields in the SCF theory. 
Considering a flow velocity v, equation (131) can be extended to the convective 
CH equation with a noise term [265]: 

^+v.V,^.V?| + <. ,132) 

Fluid dynamics is described by a continuity equation and an extended NS equation, 
which can generally be written as 

+ V-Vr^ V = -VrP + 2Vr-Avisc + Kuvi + K^t + kglast , (133) 

where p is the mass density, p is the pressure, /ivisc is the viscous stress tensor, 
and ksurf considers the surface tension from the interfaces [259]. The term kgxt 
represents external forces like gravity, while keiast ^ (^o^^(^^) ~ 1) with an ap- 
propriate length scale rg denotes the elastic response of polymers to a flow field 
leading to viscoelastic behaviour [266]. An alternative description of viscoelasticity 
is the Oldroid-B model [267], which has been used together with the PF method 
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to study the elongation and burst of viscoelastic droplets [268] and the coalescence 
of polymer drops and interfaces [261, 269]. 

The connection of the CH and the NS equation is called "model H" in the 
nomenclature of Hohenberg and Halperin [1]. In many cases, the term kgiast is 
neglected so that the polymer phase is described as a Newtonian fluid. In principle, 
all terms on the right-hand-side (besides p) can be coupled to the PF variable V^, 
since the involved material parameters may depend on the composition. Model H 
has been used to study spinodal decomposition with and without external velocity 
fields and Rayleigh- Taylor instabilities [257, 259, 270]. The described method with 
a slightly different CH free energy has been used to study the shape of fiuid films 
on a dewetting substrate [271]. 



5.2.2.3. Polymer crystals. Another important topic is the growth of polymer 
crystals in a polymer melt or solution [260, 272, 273]. So far, polymer crystals have 
not been studied explicitly with the help of the PFC method. Polymer crystalliza- 
tion has been investigated with the standard PF method, in which the extension 
of the crystalline phase is characterized by a phase field '0, while the crystalline 
structure is not considered. The dynamics of the crystal growth is described with 
the AUen-Cahn equation [274] 

where T again consists of a bulk free-energy density and a gradient term. Fre- 
quently, the heat transport is taken into account by a heat transfer equation with 
an extra term considering the production of latent heat at the crystal growth front 
[260, 273]. In various cases, two PF variables '^mat ^ind '^cryst ^ire used to repre- 
sent the concentration (or relative concentration) and the crystallization. This way 
the system is investigated with model C PF equations [258, 275]. By adding fiow 
equations, a heat equation, and external fields, like an electrical field [258], a wide 
range of multiphase polymer systems can be studied. 



5.2.2.4. Hydrogels. The PF method has also been applied to hydrogels [265], 
where the calculations consider free ions and charges fixed to the polymer network. 
The connectivity of the polymer network results in a shape elasticity as the hy- 
drogel is swollen. This model is used to study the swelling of the hydrogel as the 
ion concentration in the surrounding is changed. Its PF variable ip{r^t) changes 
continuously (and monotonously) from ip ^0 in the surrounding bath to t/^ ^ 1 in 
the hydrogel region. The dynamics of V^(r, t) is described by 

^ = -M^(h'{^)(^fisl + ^slk - (P-Pe^)) + ^g'i^P) - K^Vli!^ , (135) 
corresponding to equation (134), where is a mobility. Here, the term 

+ ^^-^4k -{P- Peq) (136) 

is the driving force for the phase field that depends on the osmotic pressure p 
and the elasticity of the hydrogel, which is given by the shear modulus and 
the difference A — A of the Lame constants inside and outside the hydrogel. The 
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deformation of the hydrogel is described by the strain 



1 f dvi dv^ 
0- 



and the osmotic pressure at equihbrium is denoted by Peq- The function h{ip) = 
i/j'^ {3 — 2t(;) in equation (135) satisfies h{0) = and h{l) — 1. Its derivative = 
dh/dtp vanishes at '0 = and ^ = 1 so that the first term in equation (135) 
contributes only at the interface between hydrogel and bath. In the second term, 
the quantities W and Vm denote the potential height and the molar volume of the 
mixture, respectively, and the function g(t(;) = '0^(1 — '0^) provides a double- well 
potential with minima at the bulk phases. Finally, the parameter controls the 
diffusion of the phase field. The sum of the second and the third term in equation 
(135) is the functional derivative of the free-energy functional in equation (128), if 
we set t/j ^ t/j - 1/2, 1/2, a W/Vm, and b 

Assuming suitably slow swelling dynamics, mechanical equilibrium provides a 
relation between strain sij and pressure p: 

^ (h{'^|,){2^ieij + \ekk5ij) + (1 - hmXeukSij - h{^){p - p'^<i)<5i,) = . (138) 

The osmotic pressure is induced by the average concentration difference of free ions 
inside and outside the hydrogel. If q is the local concentration of type i, one has 

P(i) Jdri;{r,t)ci{r,t) /dr (1 - V>(r, t))c,(r, Q 

NAkBT ^ /drV(r,i) /dr (1 - V(r, t)) ' ^ ^ 

With equations (138) and (139) the first term in equation (135) can be calculated, 
if the local concentrations q of the ions are known. The time dependence of the 
ion concentrations is given by the Nernst-Planck equation 



~dt 



= V 



(a (VrQ + l^iVA^BTc, VrC/el)) , (MO) 



where ions of type i have a concentration q, a diffusion constant Di and a valence 
Zi. The quantity is the Faraday constant and the electric potential C/gi is given 
by the Poisson equation 

Vr-(^VrC/el)+i^r( 5^^zQ + ZfCf ) = 0, (141) 



where e is the absolute permittivity of the medium and the quantities Cf and Zf 
are the concentration and the valence of the fixed charges on the polymer network, 
respectively. Since the total number of fixed charges does not change, one has 

for the concentration of fixed charges Cf , which is c^*^ at equilibrium. With the help 
of these equations, the system, consisting of a charged polymer-network in a solvent 
with ions, has been investigated and the swelling behaviour of the hydrogel has been 
analysed. This section demonstrates how the PF equation can be coupled with 
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other equations that control the mechanical, chemical, and electric properties of 
the system. It shows that short-range interactions at the moving hydrogel-solution 
interface can be coupled with long-range fields and dynamics. 



5.2.2.5. Phase-field-crystal method for polymer systems. The connectivity of 
polymer chains makes it difficult to apply the PFC method on the atomistic struc- 
ture of polymers. However, PFC models are a useful tool to investigate structure 
formation on a larger scale. One example is the investigation of density correlations 
in polymer solvents. This has been done by Praetorius and Voigt, who studied the 
density proffie of a polymer solvent in the vicinity of moving colloids [187]. The 
model considers soft interactions of polymer coils that flow around spherical col- 
loids. The colloid diameter and the effective polymer diameter are of the same 
order, leading to pronounced density oscillations in the vicinity of the colloid. The 
calculations are based on a free-energy functional of the SH-type, 

a|l = c/d.(|(-.+ (l + Vj)^). + ^+p^,), (143) 

where the phase field ^^(r, t) = 2(po — t)) characterizes the deviation of the local 
density p(r,t) from the average density po, C ^ l/(48po)5 and Ui{r,t) represents 
the rotationally-symmetric potential inferred by the colloid in units of /cbT. The 
free-energy expression is derived analogously for equation (53). The calculation of 
a polymer solution passing a fixed colloid shows pronounced density modulations 
around the colloid with pronounced oscillations in the side, where the flow hits the 
colloid. The example shows that the PFC model can also be used to analyse order 
in a fluid polymer system. 

Typically, the occurrence of structured patterns results from the interplay of 
competing interactions. In the case of crystals there are typically short-range re- 
pulsions on one hand and the external pressure plus attractive forces on the other 
hand. In various soft matter systems, other types of competing interactions may 
lead to periodic patterns on a distinctly longer scale. Seul and Andelman have 
described many pattern formation phenomena with the help of phase modulations 
with a preferred wavelength [4] . The variety of materials ranges from fluid films of 
dipolar molecules over ferrofluids to fluid membranes and block copolymers. 

We consider a system with two components A and B that tend to demix. 
The phase variable defines the deviation from the homogeneous distribution. 
A Ginzburg-Landau expansion of the free energy is given by 

jr^ = J-o[^r] + ^ y"dr (Vr^'(r))2 = To[^] + ^ jdk^'{-k)k''^ik} , (144) 

where ^(k) is the Fourier transform of ^^(r). The bulk contribution is typi- 

cally a Landau free-energy expansion. Without extra interactions, corresponds 
to the standard PF model. Now we assume a long-range interaction between 
molecules at positions r and r', which is proportional to a correlation function 
^(r,r'). This provides an extra free-energy term 

J^g^Y J^^ ydr'*(r)^(r,rO*(rO (145) 

with a constant Cg. A simple example, given in reference [4], is a 2D film of 
molecules A and B, in which molecules A have a dipole moment /i. In this case. 
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one has Cg — —ji^ and ^(r,r') = |r — r'|~^ [4]. If ^(r,r') = g{\Y — r'|), the free 
energy J^g can be expressed as 



Jdk^{-k)g{k)^{k) , (146) 



where g{k) is the 2D Fourier transform of g{\r — r'|). Combining equations (144) 
and (146), one has 



= Jo + ydk#(-k)pe(A;)*(k) (147) 



with gc{k) = Cgg{k) + hk'^. It can be shown that ^c(^) has a minimum for a finite 
fcm and that the system forms a periodic pattern with wave number k^. 

A free-energy expression hke in equation (147) is also found in the Leibler-Ohta- 
Kawasaki (LOK) theory for AB-block copolymers [276, 277]. The theory allows the 
investigation of microphase separation providing lamellar, cylindrical, and spherical 
domains. The theory starts with a model Hamiltonian, which integrates interactions 
along all pairs of polymers. The ansatz corresponds to that of the SCF theory 
[243, 278]. As a next step, the LOK theory applies a random-phase approximation 
(RPA), which leads to an expression for the free energy = ^^g + ^^hoi^j? where 
Tg corresponds to equation (146) with 

Cg5(r,r') = ^r(r,r') - 2x5{v - r') , (148) 

and J-ho[^] contains higher-order terms of ^. In equation (148), x is the Flory- 
Huggins parameter. With some approximations, a simple expression for r(k, k'), 
which is the Fourier transform of r(r,r'), is derived [277]. This leads to a free- 
energy expression of the form 

•^^=2^(2^/^'^^(-^)^^^^^^^^ ^'''^ 

with a polymerization degree N ^ the average density po? and 

h{k) = Bk^ + Ak-'^ - X . (150) 

Here, A = 3/(7V(/)^(l - and B = N/{A^{1 - </))) with (j) being the fraction of A 
monomers per polymer. Furthermore, x = 2poA^X ~ ^o(0) with a function /io(0)- 
Again, minimization of gc{k) provides a preferred wave number kc = (A/B)^/^^ 
while Xc = 2{ABY/'^ is the critical value, above which gdkc) gets negative and the 
spatially homogeneous distribution gets unstable. 
In the vicinity of the critical point, ^c(^) can be expanded around /cc, which leads 

to 



2poNkl (27r)^ 



jdk^{-k)[- klix - Xc) + B{e- klY) + ^^(k) . (151) 



Transformed into real space, one has 
B 



PoNkl 
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If we set ko = kc and /3 = k'^{x~Xc)/B, the integrand in equation (152) corresponds 
to that in the PFC free energy in equation (1), up to a term proportional to 
which is included in J-ho[*]- The free energy = J^g+Tho[^] with J^g as in equation 
(152) has already been derived by Leibler. It is only valid in the weak-segregation 
case close to the critical point. 

Tsori and Andelman have investigated a copolymer system in contact with a 
chemically patterned substrate [279]. Such systems are important for the creation 
of bulk structures controlled by a surface structure. Considering systems reasonably 
close to the critical point, the authors have used a bulk free-energy functional of the 
PFC-type. The interaction with the surface is described by a free-energy functional 



^s 



Jdr {Us{r)^{r) + Ts^\r)) . (153) 



Here, the surface potential Us{y) prefers A monomers in regions with Us{y) < 0. 
The ^^^-term of J-'^ is a surface correction to the Flory parameter x- Depending 
on whether the system is below or above the critical point the bulk phase is 
either homogeneous or lamellar. In both cases, the surface pattern modulates the 
bulk phase. The Fourier transformation of the surface pattern shows that Fourier 
modes with k > kc persist only close to the surface, while modes with k < kc can 
propagate over long distances. 

In the strong segregation regime, which is present in segregating systems far from 
the critical point, the LOK theory provides a solution in which the /c~^-term in 
equation (150) is treated more accurately, leading to 

^ = ^ /^^ (^(Vr*)' + X^' + i ydr'*(r)G(r, rO*(rO) + Jkof^] , (154) 

where the Green's function G(r, r') obeys 

-V?G(r,rO = 5(r-rO . (155) 

In the strong segregation limit the interfacial regions are small compared to the 
size of the phase domains. The LOK theory for the strong-segregation limit has 
been combined with the PF method in order to investigate fluctuations in block 
copolymer systems [280]. The dynamics of the relative concentration difference 
is studied with a CH equation as given in equation (129). As the free energy in 
equation (154) includes long-range interactions explicitly, the method goes beyond 
the standard PF theory and also beyond the standard PFC method. However, 
inserting equation (154) into equation (129) provides a rather simple equation, 
namely 

^ = V? f 2x - A-r^l) * - ^* ^ (156) 



dt ''V PoN 7 poN 

where we have used equation (155) and assumed that M is constant. 

The combination of LOK theory and PF method allows the calculation of sta- 
tistical and dynamical properties of the diblock copolymer system, down to the 
length scale of the A-B interface width. In reference [280] the method has been 
used to investigate fluctuations of the domain boundaries in a periodic structure 
of alternating A and B layers. The method can be easily transferred to other types 
of micro-segregations. 
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Another example of pattern formation in polymer system exists in the case of 
fluid membranes. Such membranes typically consist of amphiphiles, which aggre- 
gate in a layer such that the polar end groups screen the non-polar hydrocarbon 
chains from the surrounding water. In most cases, the membranes consist of two 
adjacent layers that point towards each other with their hydrophobic chains, but 
there are also various examples of fluid single- layer membranes. Seul and Andelman 
introduced a simple model of a fluctuating membrane consisting of two components 
A and B, which influence the bending behaviour of the membrane [4]. The local 
concentration is characterized by a quantity with = 1 for pure A regions 
and = for pure B regions. In the most simple case, the membrane consists of 
a monolayer. The undulations of the membrane are described by a height profile 
/i(r), where r is a 2D vector in the xi-X2-plane. For the elastic part of the free 
energy the expression 



is used. Here, a is the surface tension and is the bending modulus. The co- 
efficient A in the last term considers a coupling between the curvature and the 
local composition ^. The term Vr/i(r) is positive for convex bending and nega- 
tive for concave bending. Thus, the free energy depends on the bending direction. 
This sign-dependent bending term does not vanish if the components A and B are 
asymmetric. The total free energy of the system is given by = + J^c with 
the Ginzburg-Landau expression from equation (144). By minimizing T with 
respect to the membrane shape /i(r), the authors find the expression 



with = b — Is? The important aspect is here that the coupling of the fields 
h{v) and ^^(r) leads to a coefficient V in front of the term (Vr^)^, which may 
be negative. For V < 0, the membrane shows a curvature instability, and periodic 
patterns occur in the shape and the composition of the membrane. 
Partial integration of the free energy in equation (158) leads to 



where boundary terms have been neglected. In the relevant case of < the free 
energy can be written as 








(159) 





(160) 



with the reference wave number 




(161) 



If To is a double-well potential that can be approximated by a fourth-order poly- 
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nomial in the free energy in equation (160) corresponds to the PFC free energy 
in equation (1). As a simple example we consider that To is a symmetric potential 

•^ = «(*'-^o)', (162) 

where ^ = — is the deviation of from the overall A concentration ^. Since 
integrals over terms linear in are constant, we can write the total free energy as 

J- = 4a Jdr (|(- + R^kl,) + R^ki, + V^)')# + J*^) + (163) 



with a constant Cs and 



If we introduce the rescaled spatial coordinates r = Rr and the quantity ko = Rkj^^i^ 
the free energy is given by 

J- = AaR^ U(%(-Ul+\^\ + (kl + vi] + ^ V . (165) 



8a / V ^ V y 4 



The integral term is proportional to the free energy in equation (1) with /3 = 
^l+\h'\/M. 

5.2.3. Comparison of the methods 

The SCF theory, the PF theory and the PFC theory are three alternative methods 
for studying polymer systems. Depending on the specific problem, each of the three 
methods may be most suited. The SCF theory has originally been an equilibrium 
theory, which has been extended during the last years to study dynamic systems. 
Dynamic SCF uses a Langevin equation, in which the dynamics of the system is 
controlled by the functional derivative of a free-energy functional. The same type 
of equation, with or without a random noise term, is used in the PF and the PFC 
theories. The three methods differ in the form of the free-energy functionals. 

The SCF method is more microscopically founded than the other two. The un- 
derlying concept takes into account the behaviour of single polymer chains in an 
effective mean field provided by the surrounding polymers. This way, the connec- 
tivity and the entropy of the chains are taken into account more explicitly than 
in the PF or the PFC method. The price for this accuracy is a comparably high 
calculation effort, which limits the treatable length and time scale of the studied 
system. If details of molecular fluctuations and the connectivity of the chains are 
crucial, the SCF method is most accurate. The PF and PFC method depend on 
various parameters that have to be taken from other calculations and experiments. 
If these parameters are known, the two methods are an elegant and powerful way 
to study the dynamics of multiphase systems on a large scale. 

The SCF theory has been used to study microphase separations in copolymer 
systems. In the LOK theory, the free energy used in the SCF method is simplified 
with the help of the RPA. Close to the critical point, the LOK free-energy functional 
is of the PFC type. This way the PFC method can be used for copolymer systems 
as an approximation of the SCF theory, valid in the weak segregation limit. By 
adding a surface free-energy functional, also surface interactions have been taken 
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into account [279]. In the more usual case of strong segregation, a more general 
LOK free-energy functional has been inserted into the PF equation [280]. The 
resulting model can be seen as an extended version of the PFC method applied on 
copolymers. 

Free energies of the PFC type, which lead to pronounced periodic structures 
with a preferred wavelength, can be found for various polymer systems. As pointed 
out in reference [4], there are two mechanisms that can stabilize patterns with a 
preferred wavelength. One is the competition of long-range and short-range forces. 
The other is the interplay of two order parameters, like the shape and the com- 
position of multi-component fluid membranes. In many cases, one can derive a 
free-energy functional of the PFC type, including second- and fourth-order spatial 
derivatives of the phase field. In the case of the fluid membrane, the V^^^ term is 
proportional to the bending modulus k. Thus, by reducing the PFC equation 
can be transformed smoothly into a PF equation. At the same time, the reference 
wave number fcref diverges. This means that the equilibrium system shows a com- 
plete phase separation rather than a microphase separation. While the PFC is well 
suited for studying the dynamics of systems that form periodic structures in equi- 
librium or in a steady state, the PF method enables the study of "intermediate" 
pattern formation of phase-separating systems (where the intermediate structure 
can be extremely long- living) . 

As shown for the case of colloids flowing in a polymer solution [187], the PFC 
method can be used to consider the steric repulsion of polymer coils. Here, it reveals 
oscillatory structures with a wavelength in the range of the coil diameter, while 
similar patterns are absent in calculations comparable with the PF theory. The PF 
theory has been applied on various polymer systems, including polymer solutions, 
polymer melts, growing polymer crystals and hydrogels. It has been demonstrated 
that system properties like charge and ion distributions, flow velocity, gravitation, 
and hydrodynamics can be successfully included into the model. The according 
methods can be transferred to PFC calculations in polymer systems, resulting 
in extra terms in the free-energy functional, additional differential equations or 
advective PFC equations. 

The SCF theory is useful, if details of molecular fluctuations and connectivity 
of the chains are crucial. The PF and PFC method depend on various parameters 
that have to be taken from other calculations and experiments. If these param- 
eters are known, the methods provide an elegant and powerful way to study the 
dynamics of multiphase systems on a large scale. The application of PFC models 
on polymer systems has just started and opened many perspectives for studying 
polymer systems with pronounced pattern formation. 



5.3. Application to liquid crystals 

The different forms of PFC models for orientable particles as discussed in section 
3.1.3 can be applied to compute liquid crystalline systems under various circum- 
stances. First of all, the two-dimensional bulk calculations presented in section 
3.2.2 clearly need to be extended to three spatial dimensions. Moreover, a huge 
application area is found for confined systems. For instance, boundary conditions 
to the director fields set by external walls can lead to forced topological defects 
of the orientational order [281]. Here, PFC models for liquid crystals provide a 
flexible tool to access those numerically. Next, liquid crystalline droplets in air or 
vacuum provide another type of confinement which induces quite peculiar orienta- 
tional fields (so-called tactoids [282] or smectic droplets [282]) and again the PFC 
approach would provide a microscopic avenue to approach those. It would further 
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Figure 48. Time evolution of the Schlieren pattern oc ^^le a::i-a::2-plane, which exhibits a dynamical 

coarsening. The symbol to denotes a characteristic time. (Reproduced with permission from S. Yabunaka 
and T. Araki, Polydomain growth at isotropic-nematic transitions in liquid crystalline polymers, Phys. 
Rev. E 83 (2011), 061711, no. 6, DOT: 10. 1103/PhysRevE.83. 061711 © 2011 by the American Physical 
Society.) 

be interesting to generalize the PFC model on curved surfaces in order to explore 
the structure of thin liquid crystalline bubbles, see, e.g., references [283, 284] for 
simulation predictions. Last not least the structure of interfaces between two coex- 
isting liquid crystalline phases needs future attention, in particular for phases with 
positional order. 

A further broad range of applications have to do with dynamics. One recent ex- 
ample was performed in reference [158], where the CMA was employed to compute 
the dynamic coarsening of a disturbed nematic phase. The time evolution of such 
a nematic phase is shown in the Schlieren patterns in figure 48. 

More applications concern the dynamics of topological defects, nucleation in 
liquid crystalline systems [285], and the orientational dynamics induced by external 
switching fields [286]. Hence a flourishing future of many more PFC computations 
for liquid crystals is still lying ahead. 

6. Summary and outlook 

Motivated by the spectacular advances phase-field-crystal modelling has made in 
recent years, we have reviewed its present status. Besides presenting the original 
PFC model together with its derivation from dynamical density functional the- 
ory, we have reviewed many of its numerous extensions, including those aimed at 
describing binary solidification, vacancy transport (VPFC), anisotropic molecules 
(APFC), hquid crystals, and a quantitative description of real systems. We have 
reviewed, furthermore, a broad range of applications for metallic and soft matter 
systems (colloids, polymers, and liquid crystals), and for phenomena like the glass 
transition, and the formation of foams. We have discussed open issues such as cou- 
pling to hydrodynamics and the possibility of making quantitative PFC predictions 
for real materials. The main question at present is what further steps need yet to 
be made to turn the PFC-type models into even more potent modelling tools. 

To summarize the present state of affairs, it seems appropriate to recall some of 
the concluding remarks of a CECAM workshop dedicated to DDFT- and PFC-type 
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approaches held in 2009 in Lausanne [194]. It appears that despite the advances 
made meantime, some of the major issues identified there need yet further atten- 
dance. These are the following: 

(i) How to build numerically efficient, quantitative PFC models for a broad spec- 
trum of metallic materials ? The PFC models incorporate microscopic physics in a 
phenomenological manner. The respective local free-energy functional and the sim- 
plified dynamics lead to equations of motion that can be handled fairly efficiently 
with advanced numerical methods so that simulations containing up to a few times 
10^ particles/atoms can be performed with relative ease. A major aim here is to 
develop a methodology for tuning crystal symmetry, lattice spacing, elastic con- 
stants, surface energy, dislocation core energy, dislocation mobility, etc. without 
sacrificing numerical efficiency. Along this line, methods have been proposed for 
constructing PFC free energies that allow for the tuning of the crystal structure 
[132, 135, 136]. The amplitude equations represent an appealing alternative (see, 
e.g., references [19, 31, 131, 169]), in which the density field is expressed in terms 
of slowly-varying amplitudes, modulated by the fundamental spatial periodicity of 
particle density. As demonstrated, this approach realizes a truly multi-scale ap- 
proach to phase transitions in freezing liquids [20]. Alternatively, one can work 
directly with the scaled density field of the PFC models and introduce additional 
model parameters, which can be fitted so that a required set of physical properties 
is recovered, as done in the case of pure bcc Fe [133]. 

(a) How to construct effective, low-frequency representations from DFT/DDFT? 
Provided that one had an accurate and predictive density functional that incorpo- 
rates interaction potentials between the constituent species in a multi-component 
system, it would become possible to develop an effective description that en- 
ables quantitative simulations for microscopically-informed continuum systems 
that evolve on diffusive time scales. However, one needs to develop first such free- 
energy functionals. Next, the dynamics of the relevant degrees of freedom should 
be projected out from the full DDFT description. It may be expected on physi- 
cal grounds that the shape of a single density peak would relax much faster than 
the distance between different peak centres. Accordingly, one could "slave" the 
high-frequency modes associated with the peak shapes to the more slowly evolving 
modes with low spatial frequencies. 

(Hi) The role of fluctuations in DDFT and PFC modelling. There is a continuing 
debate about the role of noise in the DDFT- and PFC-type models [193]. Deriva- 
tions of DDFT from either the Smoluchowski level [97] or within the projection 
operator technique [74] lead to a deterministic equation of motion without any 
noise an approximation that becomes problematic near the critical point, or dur- 
ing nucleation, where the system has to leave a metastable free-energy minimum. 
In the former case fluctuations are needed to obtain the correct critical behaviour, 
whereas in the latter case, fluctuations are needed to establish an escape route of 
the system from a metastable phase. Other approaches treat fluctuations on a more 
phenomenological level. Often, however, the noise strength, though fundamentally 
correlated with the thermal energy, is treated as a phenomenological fitting param- 
eter (see, e.g., references [34, 203]). This is a fundamental problem, shared by all 
DDFT and PFC approaches. We note that the addition of noise to the equation 
of motion in continuum models is not without conceptual difficulties (see reference 
[195]), even if noise is discretised properly during numerical solution. For example, 
in the presence of noise, the equilibrium properties of the system change. Further- 
more, transformation kinetics generally depends on the spatial and temporal steps 
and in the limit of infinitely small steps in 3D the free energy of the PFC systems 
diverges, leading to an ultraviolet "catastrophe". Evidently, an appropriate "ultra- 



July 3, 2012 0:16 Advances in Physics aip 

Advances in Physics 83 

violet cutoff', i.e., filtering out the highest frequencies, is required to regularize 
the unphysical singularity. Here, a straightforward choice for the cut-off length is 
the interparticle/interatomic distance, which then removes the unphysical, small 
wavelength fluctuations [27, 28, 42, 164]. A more elegant handling of the problem is 
via renormalizing the model parameters so that with noise one recovers the "bare" 
physical properties (as outlined for the Swift-Hohenberg model in reference [197]). 
Further systematic investigations are yet needed to settle this issue. 

Appendix A. Coefficients in the PFC models for liquid crystals 

A.l. PFC model for liquid crystals in two spatial dimensions 

In the contributions (94)- (96) of the local scaled excess free-energy density, the 
coefficients 

^i=8M0(1), ^2 = -2M0(3), ^3 = ^M2(5) (A1) 

are associated with a gradient expansion of ^^(r). These coefficients also appear in 
a different form in the original PFC model [26]. The further coefficients are given 
by [147, 287] 

Si = 4(Mii(2)-M?(2)), (A2) 
B2 = 2{M\{2)-Ml^{2)), (A3) 
S3 = -M2 2(3)-M0(3), (A4) 



Ci=4MJ(1), C2 = mJ(3)-^M^2(3), C3 = -ML2(3), (A5) 



and 



Di = 2M2(1), D2 = -MI{3). (A6) 



So far, all these coefficients can also be obtained by using the second-order Rama- 
krishnan-Yussouff functional for the excess free energy. The remaining coefficients, 
however, result from higher-order contributions in the functional Taylor expansion 
[147]. In third order, one obtains for the homogeneous terms the coefficients 

El = 32^41, (A7) 

E2 = 16(Moo'' + 2M[]i), (AS) 

£;3 = 8(Moo' + 2M02), (A9) 

£;4 = 8(2Moo21 + mJJ) (AIO) 
and for the terms containing a gradient the coefficients 

Fi = 16(m^I'' - 2M0^i + Mg?) , (All) 
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F2 = 16 (Mq^i - M^r' + Mgi - Mjr') , (A12) 
Fs = -16 (Mq 20 _ - Mgi + Mj?) , (A13) 

F4 = -8(Moi'"' - 2 Mo/' + Mjr') , (A14) 
F, = -4(Mo-2-i - - + m2-i) , (A15) 
Fe = 8 (Mof - + _ ^2-2 j _ (^^g) 

In fourth order, only homogeneous terms are kept. The corresponding coefficients 
are 

G, = 128MZ, (A17) 
G2 = 192(m^,\^' + mZ), (A18) 
G3 = 96(m^^^ + mZ), (A19) 
G4 = 96 (2 Moo^oi + Moo^ii + mOJi) , (A20) 

G5 = 48 (Moo'o'' + Mooo ') > (A21) 
Gq = 48M^^\'\ (A22) 
G7 = 12Moo222_ (^23) 

All the coefficients from above are linear combinations of moments of the Fourier 
expansion coefficients of the direct correlation functions. These moments are de- 
fined through 

(«") = TT^^+Vrel' ( n /dr. rf^ cl:+l)(r") (A24) 



with the multi-index notation = (Xi, . . . , X^) for X G {/, m, r, 1, a, 0, 0r} and 
the abbreviations Mp?" = Mp?"(l^) and M^^/^^ = M[^/^^^(l, 2). The expansion 
coefficients of the direct correlation functions are given by 



- ld<P^ c("+')(r", 0")e-('"-'^^+™°-<^") (A25) 

with ri - r^+i = i?^u((/PR.), = u(^z), = ^1 - ^r,, and 0^ = (^1 - (pi^i. 

When the system is apolar, the modulus P(r) of the polarization P(r) is zero 
and its orientation p(r) is not defined, while the direction n(r) associated with 
quadrupolar order still exists. Then, symmetry considerations lead to the equalities 

c%{R) = cg(i?) , c%2{R) = cg(^) > - 4?(^) (A26) 



between expansion coefficients of the direct pair-correlation function and to the 
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equations 

M^i(2) = M?(2) , M^i(2) = M{(2) , M^2(2) = M^(2) (A27) 

for the generalized moments. A consequence of these equations is that the coeffi- 
cients Bi and B2 vanish and becomes more simple. 



A. 2. PFC model for liquid crystals in three spatial dimensions 

The coefficients in equation (105) appear in three different groups. The first group 
consists of the three coefficients 

Ai = 8 0ooo(0) , ^2 = - ^ f^ooo(2) , ^3 = ^ f^ooo(4) , (A28) 

that are already known from the original PFC model [26] and belong to the gradient 
expansion of the translational density. In the next group, the two coefficients 

Bi = ^ f^220(0) , ^2 = ^^022(2) , (A29) 

that go along with the nematic tensor and the coupling of its gradient with the 
gradient of the translational density, are collected. The last group contains the 
Frank constants 

that appear in the Frank free-energy density [145, 157]. All these coefficients are 
expressed in terms of the generalized moments 

POO 

^hiAn) = TT^/'Pr'ef / dv T^^^uiMir) (A31) 

^0 

with the expansion coefficients 

I — ^— — min{/i,/2} 

^hiAr) = y^j^ C'(/i,/2,/,m,-m,o) 

m=-min{/i,/2} (A32) 

X ydui Jdu2 c(^)(re3,ui,U2)Y/^^(ui)Y/^_^(u2) , 

where the symbol C(Zi, Z25 Z, mi, m2, m) denotes a Clebsch-Gordan coefficient [288], 
Y/^(u) is a spherical harmonic, 63 stands for the Cartesian unit vector co- 
directional with the X3-axis, and ~ denotes complex conjugation. 

As before, equalities between generalized moments with different index- 
combinations, that can be derived using symmetry considerations [59], have been 
taken into account in order to reduce the set of generalized moments in equation 
(105) to its seven independent members. 
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2D 


two spatial dimensions 


OTA 

3D 


three spatial dimensions 


1M-FJ:^C model 


Single- mode Pl^ C model 


2M-PFC model 


two-mode PFC model 


APJ:^C model 


anisotropic Pf^ G model 


AiG instability 
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Asaro- liller-Grmield instability 


bcc crystal structure 


body-centred cubic crystal structure 


bet crystal structure 


body-centred tetragonal crystal structure 


"DA /"ID 


boundary value problem 


CH equation 


Cahn-Hilliard equation 




XX 1 •! 'x • X • 

constant-mobility approximation 


TATA 77 
DDr^ 1 


dynamical DFT 


DbP method 


1 • X 1 xx'l xll 

dynamic external potential method 
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density functional theory 


DLVU potential 


Derj aguin-Landau- Verwey-Overbeek potential 


DMD simulation 


diffusive MD simulation 


DSCF theory 


dynamic SCF theory 


EAP-MD simulation 


embedded-atom-potential MD simulation 


ELE 


Euler-Lagrange equation 


EOF-PFC model 


eighth-order fitting PFC model 


EOM 


equation of motion 


fee crystal structure 


face-centred cubic crystal structure 


FD scheme 


finite-difference scheme 


FMT 


fundamental-measure theory 


GRP-PFC model 


Greenwood-Rottler-Provatas PFC model 


hep crystal structure 


hexagonal close packed crystal structure 


HS potential 


hard-sphere potential 


LJ potential 


Lennard- Jones potential 


LOK theory 


Leibler-Ohta-Kawasaki theory 


MC simulation 


Monte-Carlo simulation 


MCT 


mode-coupling theory 


MD simulation 


molecular-dynamics simulation 


Model A 


relaxational dynamical equation for a non-conserved order 




parameter 


Model B 


relaxational dynamical equation for a conserved order 
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Model C 

Model H 
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NS equation 
Oldroid-B model 
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RPA 
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and a conserved order parameter 

dynamical equation for a conserved order parameter coupled 
to flow 

modified PFC model 
Navier-Stokes equation 

constitutive model for the flow of viscoelastic fluids 
phase-field-crystal model 

dynamical equation for the original PFC model without CMA 

dynamical equation for the original PFC model with CMA 

phase-field model 

reciprocal lattice vector 

random-phase approximation 

simple cubic crystal structure 

self-consistent field theory 

Swift-Hohenberg model 

vacancy PFC model 



